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When ambient seismic waves pass near and under an interferometric gravitational-wave detector, 
they induce density perturbations in the earth, which in turn produce fluctuating gravitational forces 
on the interferometer's test masses. These forces mimic a stochastic background of gravitational 
waves and thus constitute a noise source. This seismic gravity-gradient noise has been estimated and 
discussed previously by Saulson using a simple model of the earth's ambient seismic motions. In this 
paper, we develop a more sophisticated model of these motions, based on the theory of multimode 
Rayleigh and Love waves propagating in a multilayer medium that approximates the geological 
strata at the LIGO sites (Tables 0-|l\j), and we use this model to revisit seismic gravity gradients. 
We characterize the seismic gravity-gradient noise by a transfer function, T{f) = x{f) /W (f), from 
the spectrum of rms seismic displacements averaged over vertical and horizontal directions, W{f), 
to the spectrum of interferometric test-mass motions, x{f) = Lh{f); here L is the interferometer 
arm length, h{f) is the gravitational-wave noise spectrum, and / is frequency. Our model predicts 
a transfer function with essentially the same functional form as that derived by Saulson, T ~ 
47rGp(27r/)~^/3(/), where p is the density of the earth near the test masses, G is Newton's constant, 
and /3(/) = 7(/)r(/)/3'(/) is a dimensionless reduced transfer function whose components 7 ~ 1 
and r ~ 1 account for a weak correlation between the interferometer's two test masses (Fig. |l|) and 
a slight reduction of the noise due to the height of the test masses above the earth's surface. This 
paper's primary foci are (i) a study of how f3'{f) ~ /3(/) depends on the various Rayleigh and Love 
modes that are present in the seismic spectrum (Figs. and Table ^, (ii) an attempt to estimate 
which modes are actually present at the two LIGO sites at quiet times and at noisy times, and (iii) 
a corresponding estimate of the magnitude of /?'(/) at quiet and noisy times. We conclude that at 
quiet times ~ 0.35-0.6 at the LIGO sites, and at noisy times j3' ~ 0.15-1.4. (For comparison, 
Saulson's simple model gave (5 — P' — l/\/3 = 0.58 .) By folding our resulting transfer function into 



the "standard LIGO seismic spectrum" [Eq. (1.29)], which approximates W{f) at typical times, we 



obtain the gravity-gradient noise spectra shown in Fig. 2. At quiet times this noise is below the 
benchmark noise level of "advanced LIGO interferometers" at all frequencies (though not by much at 
~ 10 Hz); at noisy times it may significantly exceed the advanced noise level near 10 Hz. The lower 
edge of our quiet-time noise constitutes a limit, beyond which there would be little gain from further 
improvements in vibration isolation and thermal noise, unless one can also reduce the seismic gravity 
gradient noise. Two methods of such reduction are briefiy discussed: monitoring the earth's density 
perturbations near each test mass, computing the gravitational forces they produce, and correcting 
the data for those forces; and constructing narrow moats around the interferometers' corner and end 
stations to shield out the fundamental-mode Rayleigh waves, which we suspect dominate at quiet 
times. 
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I. INTRODUCTION AND SUMMARY 

Now that the LIGO/VIRGO international network of 
gravitational-wave detectors is under construction, 
it is important to revisit the various noise sources that 
will constrain the network's ultimate performance. Im- 
proved estimates of the ultimate noise spectra are a foun- 
dation for long-term planning on a number of aspects 
of gravitational-wave research, including facilities design, 
interferometer R&D, data analysis algorithm develop- 
ment, and astrophysical source studies. 

In this paper and a subsequent one |^ we revisit 
gravity- gradient noise — noise due to fluctuating Newto- 
nian gravitational forces that induce motions in the test 
masses of an interferometric gravitational-wave detector. 
Gravity gradients are potentially important at the low 



end of the interferometers' frequency range, / < 20 Hz. 
Another noise source that is important at these frequen- 
cies is vibrational seismic noise, in which the ground's 
ambient motions, filtered through the detector's vibra- 
tion isolation system, produce motions of the test masses. 
It should be possible and practical to isolate the test 
masses from these seismic vibrations down to frequencies 
as low as / ~ 3 Hz but it does not look practical to 
achieve large amounts of isolation from the fluctuating 
gravity gradients. Thus, gravity gradients constitute an 
ultimate low-frequency noise source; seismic vibrations 
do not. 

Gravity gradients were first identified as a potential 
noise source in interferometric gravitational-wave detec- 
tors by Rai Weiss in 1972 0. The first quantitative anal- 
yses of such gravity-gradient noise were performed by Pe- 
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ter Saulson [g[ and Robert Spero in the early 1980s. 
There has been little further study of gravity-gradient 
noise since then, except for some updating in Saulson's 
recent monograph |10[ . 

In his updating, Saulson concluded that the most se- 
rious source of gravity-gradient noise will be the fluc- 
tuating density of the earth beneath and near each of 
the interferometer's test masses. These density fluctu- 
ations are induced by ambient seismic waves that are 
always present; their resulting gravitational forces are 
called seismic gravity-gradient noise. Saulson p|jlO|] also 
estimated the gravity gradient noise from atmospheric 
fluctuations, concluding that it is probably weaker than 
that from earth motions. Spero |^ showed that gravity- 
gradient noise due to jerky human activity (and that of 
dogs, cattle, and other moving bodies) can be more seri- 
ous than seismic gravity-gradient noise if such bodies are 
not kept at an adequate distance from the test masses. 
We shall revisit seismic gravity-gradient noise in this pa- 
per, and gravity gradients due to human activity in a 
subsequent one ||^; Teviet Creighton at Caltech has re- 
cently initiated a careful revisit of gravity gradient noise 
due to atmospheric fluctuations. 

Our detailed analysis in this paper reveals a level 
of seismic gravity-gradient noise that agrees remarkably 
well with Saulson's much cruder estimates. Our analy- 
sis reveals the uncertainties in the gravity gradient noise, 
the range in which the noise may vary from seismically 
quiet times to noisy times, the dependence of the noise 
on the various seismic modes that are excited, and the 
characteristics of the modes that the geological strata 
at Hanford and Livingston are likely to support. The 
dependence of the noise on the modes and the character- 
istics of the expected modes are potential foundations for 
methods of mitigating the seismic gravity gradient noise, 
discussed in our concluding section. 

A preliminary version of this paper [prf was circu- 
lated to the gravitational-wave-detection community in 
1996. That version considered only fundamental-mode 
Rayleigh waves (which we suspect are responsible for the 
dominant seismic gravity-gradient noise at quiet times), 
and (as Ken Libbrecht pointed out to us) it contained a 
serious error: the omission of the "surface-source" term 
[denoted in Eq. (1.24) below] for the gravity-gradient 



force. It also contained errors in its two-geological-layer 
analysis for the LIGO Hanford site. These errors have 
been corrected in this final version of the manuscript, 
and the analysis has been extended to include more re- 
alistic models of the geological strata at the two LIGO 
sites and to include higher-order seismic modes. 

As we were completing this manuscript, we learned of 
a paper in press ||l^ by Giancarlo Cella, Elena Cuoco and 
their VIRGO-Project collaborators, which also analyzes 
seismic gravity-gradient noise in interferometric gravita- 
tional wave detectors. That paper is complementary to 
ours. Both papers analyze the RF mode (which we sus- 
pect is the dominant contributor to the seismic gravity- 
gradient noise at quiet times), obtaining the same re- 



sults in the 3-30 Hz band when the effects of geological 
stratification are neglected. But, whereas our paper car- 
ries out an extensive study of the effects of stratification 
and other modes, the Cella-Cuoco paper extends the un- 
stratified RF-mode analysis to frequencies below 3 Hz 
and above 30 Hz, and computes (and finds to be small) 
the gravity gradient noise caused by seismically-induced 
motions of the experimental apparatus and its massive 
physical infrastructure in the vicinity of the VIRGO test 
masses. 



This paper is organized as follows: In Sec. I A, we de- 



scribe the phenomenology of the seismic- wave modes that 
can contribute to ambient earth motions at horizontally 
stratifi ed s ites like LIGO-Hanford and LIGO-Livingston. 
In Sec. IB, we introduce the transfer function T(/) used 
to characterize seismic gravity-gradient noise, we break 
it down into its components [most especially the reduced 
transfer function /?'(/)], and we express it as an inco- 
herent sum over contributions from the various seismic 
modes. In Sec. we briefly describe Saulson's com- 
putation of the reduced transfer function, and then in 
Sec. |[ D| we describe ou r own computation and results. 
More specifically, in I D we gather together and summa- 
rize from the body of the paper our principal conclusions 
about f3' for the various modes at the two LIGO sites, 
we discuss the evidence as to which modes actually con- 
tribute to the noise at quiet times and at noisy times, and 
we therefrom estimate the net values of at quiet and 
noisy times. We then fold those estimates into the stan- 
dard LIGO seismic spectrum to get spectral estimates of 
the seismic gravity-gradient noise (Fig. 2). 

The remainder of the paper (summarized just before 
the beginning of Sec. |l^ presents our detailed models 
for the geological strata at the two LIGO sites, and our 
analyses of the various seismic modes that those strata 
can support and of the seismic gravity-gradient noise pro- 
duced by each of those modes. 



A. Phenomenology of ambient seismic motions in 
the LIGO frequency band 

Seismic motions are conventionally decomposed into 
two components pl|"p^: P-waves and S-waves. P-waves 
have material displacements along the propagation direc- 
tion, a restoring force due to longitudinal stress (pressure 
— hence the name P-waves) , and a propagation speed de- 
termined by the material's density p and bulk and shear 
moduli K and fi: 



cp = 



'X + 4^/3 



(1.1) 



S-waves have transverse displacements, restoring force 
due to shear stress, and propagation speed 



cs 




(1.2) 
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Here v is the material's Poisson ratio 



(1.3) 



Near the earth's surface, where seismic gravity-gradient 
noise is generated, these speeds are in the range cp ~ 
500-2000 m/s and eg ~ 250-700 m/s. However, some of 
the modes that may contribute to the noise extend down 
to much greater depths, even into the bedrock where 
CP - 5000-6000 m/s and cs ~ 3200 m/s. 

The ambient seismic motions are a mixture of P- 
waves and S- waves that propagate horizontally ( "surface 
waves"), confined near the earth's surface by horizontal 
geological strata. Depending on the mode type and fre- 
quency, the horizontal propagation speed cu can range 
from the surface layers' lowest S-speed to the bedrock's 
highest P-spccd: 250 m/s < ch < 6000 m/s. 

P- and S-waves are coupled by geological inhomo- 
geneities (typically discontinuities at geological strata) 
and by a boundary condition at the earth's surface. 
At both LIGO sites the strata are alluvial deposits 
above bedrock, with discontinuities that are horizontal 
to within 2 degrees (more typically to within less than 
1 degree). Throughout this paper we shall approximate 
the material as precisely horizontally stratified. 

Seismic gravity-gradient noise is a potentially serious 
issue in the frequency band from / ^ 3 Hz (the low- 
est frequency at which mechanical seismic isolation looks 
practical) to / ~ 30 Hz; c/. Fig. 2 below. In this fre- 
quency band, the wavelengths of P- and S-waves are 

(cp/1000ms~i) (c5/500ms-i) 
Ap = 100m , p i^r.^^ \ : As = 50m- 



(//lOHz) 



(//lOHz) 



(1.4) 



Neglecting coupling, the amplitudes of these waves at- 
tenuate as exp(— Trr/QA), where r is the distance the 
waves have propagated and Q is the waves' quality factor. 
The dominant dissipation is produced by the waves' shear 
motions and can be thought of as arising from an im agi- 
nary part of the shear modulus in expressions (IT) and 
( |l.2D for the propagation speeds cs and cp (and thence 
also from an imaginary part of the propagation speeds 
themselves). Since the restoring force for S-waves is en- 
tirely due to shear, and for P-waves only about half due 
to shear, the S-waves attenuate about twice as strongly 
as the P-waves. The measured Q-factors for near-surface 
materials are Qs 10-25, Qp - 20-50 |2|j2|l, corre- 
sponding to amplitude attenuation lengths 



Hi 



QsXs 



1000 m 



250 m 



(Qp/30)(cp/1000ms-i) 
//lOHz 
(Q5/15)(cs/500ms-i) 
//lOHz 



(1.5) 



For bedrock (and basalt that overlies it at Hanford), the 
Q's and attenuation lengths can be higher than this — 
Qp as high as a few hundred P7| . 



Shallowly seated wave modes which cause ambient seis- 
mic motions in our band, i.e., modes that are confined to 
the alluvia so ch 2500 m/s (and more typically ^ 1000 
m/s), must be generated in the vicinity of the interfer- 
ometers' corner and end stations by surface sources such 
as wind, rain, and human activities (automobile traf- 
fic, sound waves from airplanes, etc.); their attenuation 
lengths are too short to be generated from further than 
a kilometer or so. Deep seated modes that reach into the 
bedrock could originate from rather further away — at 
10 Hz and in a layer that has Qp ~ 100, cp ~ 5500 m/s, 
modes can propagate as far as ~ 20 km. 

In horizontally stratified material, the wave compo- 
nents that make up each mode all propagate with the 
same angular frequency uj = 27r/, horizontal wave vec- 
tor k ~ kk (where k is their horizontal direction, and 
k = 2tt / X their horizontal wavenumber), and horizontal 
phase speed ch = to/k. Their vertical motions differ from 
one horizontal layer to another and from P-componcnt to 
S-component. The horizontal dispersion relation uj^k) [or 
equivalently c/f (/)] depends on the mode (Figs. ||, y, ||, 
and |l^ below) . 

Geophysicists divide these surface normal modes into 
two types ||,||: 

• Love modes, which we shall denote by L. These 
are S-waves with horizontal displacements ("SH- 
waves") that resonate in the near-surface strata. 
They involve no P-waves and thus have no com- 
pression and no density variations; therefore, they 
produce no fluctuating gravitational fields and no 
seismic gravity-gradient noise. 

• Rayleigh modes, which we shall denote by R. These 
are combinations of S-waves with vertical displace- 
ments and P-waves ("SV- waves") that are cou- 
pled by the horizontal discontinuities at strata in- 
terfaces, including the earth's surface. Rayleigh 
modes are the producers of seismic gravity-gradient 
noise. 

We shall divide the Rayleigh modes into two groups: 
the fundamental Rayleigh mode, which we denote RF, 
and Rayleigh overtones (all the other modes). Rayleigh 
overtones require stratification of the geological structure 
in order to be present; they essentially consist of coupled 
SV- and P-waves which bounce and resonate between 
the earth's surface and the interfaces between strata. We 
shall further divide the Rayleigh overtones into two broad 
classes: those that are composed predominantly of SV- 
waves, denoted RS, and those composed predominantly 
of P-waves, denoted RP. In the geophysics literature, the 
modes we identify as RP are sometimes referred to sim- 
ply as P-modes, and our RS modes are referred to as 
the Rayleigh overtones. However, when RP modes are 
intermixed with RS modes in the {ch, f) space of disper- 
sion relations (as turns out to be the case at Hanford; 
cf. Fig. ^ below), a given Rayleigh overtone will continu- 
ously change character from RS to RP. Because this will 
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be quite important for the details of the seismic gravity- 
gradient noise, we prefer to emphasize the similarities of 
the two mode types by designating them both as Rayleigh 
overtones and denoting them RS and RP. 

We shall append to each Rayleigh overtone an integer 
that identifies its order in increasing horizontal speed ch 
at fixed frequency /. Each successive Rayleigh mode, 
RF, RSI, RS2, . . . (and, as a separate series, RF, RPl, 
RP2, . . .) penetrates more deeply into the earth than the 
previous one. In our frequency band, the fundamental 
RF is typically confined to within ^ As/tt ^ 10 m of the 
earth's surface. 

The RF mode is evanescent in all layers (except, at low 
frequencies, in the top layer). The overtones RSI, RS2, 
. . . are composed primarily of SV-waves that propagate 
downward from the earth's surface, reflect off some in- 
terface, return to the surface and reflect back downward 
in phase with the original downward propagating waves, 
thereby guaranteeing resonance. On each reflection and 
at each interface between layers, these modes generate 
a non-negligible admixture of P-waves. The RP over- 
tones are similar to RS, but with the propagating and 
reflecting waves being largely P with some non-negligible 
accompanying SV. 

Dissipation will cause an overtone's waves to damp out 
with depth. If that damping is substantial in traveling 
from the surface to the reflection point, the overtone will 
not resonate and will be hard to excite. Roughly speak- 
ing, the amount of amplitude decay in traveling from the 
surface to the reflection point and back to the surface is 
utt/Q where n is the mode number (or equiv alen tly the 
number of round-trip wavelengths); c/. Eqs. (1.5). The 
round-trip damping therefore exceeds 1/e for mode num- 
bers n ^ Qs/t^ ~ 5 for RS modes and n ^ Qp/tt ~ 10 
for RP modes. Correspondingly, in this paper we shall 
confine attention to modes with mode numbers n ^ 10. 

The RP modes are harder to analyze with our formal- 
ism than RS modes — typically, when RP modes turn 
on, there are many modes very closely spaced together 
and it is difficult to distinguish them. For this reason, we 
shall study only the lowest one at each site, RPl, plus 
RP modes that travel nearly horizontally in the several- 
km thick basalt layer at Hanford. We expect RPl to be 
typical of other low-order RP modes, and the basalt-layer 
RP waves to be typical also of such waves propagating 
nearly horizontally in the bedrock. 



B. Transfer functions and anisotropy ratio 

Following Saulson ||l^ , we shall embody the results of 
our gravity-gradient analysis in a transfer Junction 



Tif) ^ 



W{f) 



(1.6) 



from seismic-induced earth motions W{f) to differential 
test-mass motion i(/). The precise definitions of W{f) 



and i(/) are as follows: 

We shall denote the square root of the spectral density 
(the "spectrum") of the earth's horizontal surface dis- 
placements along some arbitrary horizontal direction by 
X{f) (units m/v'Hz), where / is frequency. We assume 
that X{f ) is independent of the chosen direction, i.e. the 
seismic motions are horizontally isotropic. This is justi- 
fied by seismometer measurements at the LIGO sites be- 
fore construction began |l3|Jl^ ] and by rough estimates of 
the diffractive infiuence of the constructed facilities (Sec. 
^). We shall denote the spectrum of vertical displace- 
ments at the earth's surface by Z{f). The quantity W[f) 
that appears in the transfer function is the displacement 
rms-averaged over 3-dimensional directions: 



W{f) = 



(1.7) 



The other quantity, i(/), which appears in the trans- 
fer function (1.6) is related to the interferometer's 
gravitational-wave strain noise spectrum h{f) by x{f) = 
h{f)L, where L is the interferometer arm length (4 km for 
LIGO). Physically, x{f) is the spectrum of the interfer- 
ometer's arm-length difference and is called the interfer- 
ometer's "displacement noise spectrum". Since x{f) and 
W{f) both have units of m/\/Hz, the transfer function 
T(/) is dimensionless. 

In this paper we shall express T(f) in terms of a dimen- 
sionless correction /?(/) to a simple and elegant formula 
that Saulson M derived: 



T(/) ^ ^/^ 



AirGp 



W{f) 



AnGp 



-Af) 



2/3(/) at/>3Hz. 



(1.8) 



Here p ~ 1.8 g/ cm is the mass density of the earth in the 
vicinity of the interferometer, G is Newton's gravitational 
constant, w = 27r/ is the angular frequency of the seismic 
waves and their fluctuating gravitational forces, and wq ~ 
2tt rad/s and r ~ 10^ s are the angular frequency and 
damping time of the test mass's pendular swing. We 
shall call /9(/) the reduced transfer function. Saulson's 
estimate for (3{f) was 



/5sa 



ilsoii 



l/\/3 = 0.58 



(1.9) 



cf. Eq. (21) of Ref. Q. Our analyses (below) suggest 
that at quiet times (3 may be ~ 0.35 to 0.6, and at noisy 
times, (3 ~ 0.15 to 1.4. Thus, Saulson's rough estimate 
was remarkably good. 

Each mode of the earth's motion will contribute to the 
transfer function, and since the relative phases of the 
modes should be uncorrelated, they will contribute to 
/3(/) in quadrature: 



/3 



(1.10) 
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The sum runs over all Rayleigh and Love modes, J € 
(RF, RSn, RPn, Ln); (3j{f) is the reduced transfer func- 
tion for mode J, with 



hn = 



(1.11) 



because the Love modes produce no gravity-gradient 
noise. The weighting factor wj is the fractional contribu- 
tion of mode J to the mean square seismic displacement 
W^, and correspondingly the wj's are normalized by 



1 



(1.12) 



Besides this normalization condition, there is another 
constraint on the weighting factors wj: each mode (at 
each frequency) has its own ratio Aj of vertical to hori- 
zontal displacement at the earth's surface: 



AAf) 



Zjjf) 
Xj{f) 



(1.13) 



We shall call this ratio the mode's anisotropy ratic^. 
Since the Love modes have purely horizontal motions, 
their anisotropy ratios vanish: 



At 







(1.14) 



It is straightforward to show that the anisotropy ratios 
for the various modes combine to produce the following 
net anisotropy in the earth's surface displacement: 



A: 



(1.15) 



At quiet times, measurements show this to be near unity 
at Hanford Q, and ~ 0.6 at Livingston while at 
noisy times it can fluctuate from ~ 0.2 to ~ 5. The mea- 
sured value of this ratio is an important constraint on 
the mixture of modes that produces the observed seis- 
mic noise and thence on the net reduced transfer func- 
tion. For example, if the observed noise is due to one 
specific Rayleigh mode J with large anisotropy ratio Aj, 
accompanied by enough Love waves to reduce the net 
anisotropy ratio to Anct — 1.0 (Hanford) or 0.6 (Liv- 
ingston), then Eqs. ( 1.10D -( i.l5| ) imply that the net re- 
duced transfer function for the seismic gravity gradient 
noise is 



/3jL =/3j1 



' 1 + 2/A^j 
1 + 



(1.16) 



In Appendix ^ it is shown that for each mode J, the 
reduced transfer function f3j can be split into the product 
of three terms: 



*Geophysicists use the name spectral ratio for 1/A 
l/(anisotropy ratio). 
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FIG. 1. The function 7(1/) that accounts for correla- 
tions of seismic gravity-gradient noise in the two corn er tes t 
masses. This function is given analytically by Eq. (L.19), 
and it appears in all of the reduced transfer functions: 

l3jif)=P'AfM27vfl/cHj)TM). 



(1.17) 



The first term, jj, accounts for the correlation between 
the gravity-gradient noise at the interferometer's two cor- 
ner test masses. It is a universal, mode- independent func- 
tion of the waves' horizontal phase shift in traveling from 
one test mass to the other: 



7j = ^{ujI/chj) 



(1.18) 



Here oj = 27r/ is the waves' angular frequency, ^ ^ 5 m is 
the distance between the two corner test masses, ch.j is 
the horizontal phase speed cr for mode J, and uj / ch = k 
is the mode's horizontal wave number. For frequencies 
and modes of interest to us, the argument y — ujI/chj of 
7 is of order unity. The function 7(2/), given by 



27r 



, , / 1 /"" , , / COS( 

7(y) = 'Wl + ^y cos0sm0cos ( j/- 



sm ( 



V2 



d4> , 



(1.19) 



is plotted in Fig. As Fig. |l| shows, 7(1/) is within about 
10 per cent of unity for all frequencies, so we shall regard 
it as unity througout the rest of this manuscript, except 
in Appendix |A[ 

The second term, Fj, in Eq. (1.17) for (ij describes 



the attenuation of the gravity gradient noise due to the 
height Ti. of the test masses above the earth's surface. We 
show in Appendix ^ that 



Fj = exp(-Lj7i/c/fj) 



(1.20) 
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For LIGO interferometers H is about 1.5 m, the fre- 
quency of greatest concern is / = io/2-K ~ 10 Hz (c/. 
Fig. 2 below), and at quiet times the dominant contri- 
bution t o the noise probably comes from the RF mode 
(c/. Sec. |[D j ) for which, near 10 Hz, ch — 330 m/s (c/. 
Figs. ^ and ||); correspondingly, Trf — 0.75. For other 
modes, ch will be larger so Fj will be closer to unity 
than this. For this reason, throughout the rest of this 
paper, except in Appendix we shall approximate Fj 
by unity. With 7j and Fj both approximated as unity, 
we henceforth shall blur the distinc tion b etween (3j and 
(3'j, treating them as equal [cf. Eq. ( 1.17 )]. 

In Appendix ^ we derive expressions for the reduced 
transfer function /3j(/) and the anisotropy ratio Aj in 
terms of properties of the eigenfunctions for mode J: de- 
note by ^Hj and ^yj the mode's complex amplitudes at 
the earth's surface {z = 0) for horizontal displacement 
and upward vertical displacement, so the mode's surface 
displacement eigenfunction is 



i{k-x—ujt) 



(1.21) 



where is the unit vector pointing downward and k — 
k/k is the unit vector along the propagation direction. 
Also, denote by Tlj{z) the mode's amplitude for the frac- 
tional perturbation of density Sp/ p at depth z below the 
surface, so 



P 



(1.22) 



Here the term ^v.J^i^) accounts for the mass moved 
above z = by the upward vertical displacement ^y. 
Then, we show in Appendix |^ [Eq. (A?)] that 



= i = ££:I = lfi (cg,//1000ms-i 
k 2^/ ^™ (//lOHz) 



(1.25) 



This has a simple explanation: (i) to produce much grav- 
itational force on a test mass, a compressed bit of matter 
must reside at an angle a > 7r/4 to the vertical as seen 
by the test mass, and (ii) bits of matter all at the same 
a 7r/4 and at fixed time have fractional compressions 



Sp/p that oscillate with depth z as e*' 



^ikz tan a 



and 



that therefore tend to cancel each other out below a depth 
l/{k tan a) ^ 1/fc . 

From Eq. (|l.24| ) we can estimate the magnitude of the 
reduced transfer function. The mode's fractional density 
perturbation TZj is equal to the divergence of its displace- 
ment eigenfuction (aside from sign), which is roughly 
k£.Hj and often docs not vary substantially over the shal- 
low depths z ^ 2sgg where the gravity gradients orig- 



inate. Correspondingly, the integral in Eq. ( 1.24 ) is ^ 

CffJ, so ^ V^MiHJ + ^Vj\yi\iHj\' + \^Vj\^) - 1, 

since the horizontal and vertical displacements are com- 
parable. 

As we shall see in Sees. H and |III B 3| below, for RP 
modes the gravity gradients produced by the surface 
and subsurface sources tend to cancel (a consequence of 
mass conservation), so /3' actually tends to be somewhat 
smaller than unity. 



/3rp ^ 0.15 



(1.26a) 



while for RF and RS modes, the surface source tends to 
dominate, so 



1 



RF 



V2 V 1 



2^2 



^.0.7, 



(1.26b) 



Aj 



V2. 



VJ\ 



HJ\ 



(1.23) 



where the ^/2 comes from the fact that when this mode is 
incoherently excited over all horizontal directions fc, its 
rms horizontal amplitude along any chosen direct ion is 
\Chj\/V^- Similarly, we show in Appendix |^ [Eq. ( |A2l| )] 
that 



3/2 



\^HJ\' + \^VJ\ 



Cvj+ / njiz)e-'''-dz 



(1.24) 



where k = lo/ch j for mode J. We shall refer to the 



^yj term in Eq. ( 1.24 ) as the surface source of gravity 
gradients, and the jTZje^'^'^dz term as the subsurface 
source. 

Note that the influence of a given density perturbation 
dies out as e~*'^, so unless TZj{z) increases significantly 
with depth, the seismic gravity gradients arise largely 
from depths shallower than the gravity- gradient e-folding 
length 



If we had normalized our transfer function to the vertical 
displacement spectrum |.^(/)| instead of the direction- 



averaged spectrum |VF(/)| [Eq. ([l.6[)], then for modes in 
which the surface source strongly dominates, /3j would 
be l/-\/2 independently of the mode's anisotropy ratio. 
In Sees. [II and IV and associated Appendices, we 



shall derive, for each low-order Rayleigh mode at Han- 
ford and Livingston, the reduced tran sfer function (3'j and 
the anisotropy ratio Aj. In Sec. ID, we shall discuss the 
likely and the allowed weight ings wj of t he va rious modes 
[subject to the constraints (|1.12D and (1.15)], and shall 



estimate the resulting net reduced transfer functions /?(/) 
for the two sites and for quiet and noisy times (Table |). 

Henceforth wc typically shall omit the subscript J that 
denotes the mode name, except where it is needed for 
clarity. 



C. Saulson's analysis and transfer function 

In his original 1983 analysis of seismic gravity-gradient 
noise [|^, Saulson was only seeking a first rough estimate, 
so he used a fairly crude model. He divided the earth 
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near a test mass into regions with size Ap/2 (where Ap is 
the wavelength of a seismic P-wave) , and he idealized the 
masses of these regions as fluctuating randomly and inde- 
pendently of each other due to an isotropic distribution 
of passing P-waves. Saulson's final analytic result [his 
Eq. (21)] was the transfer function (1.6) with j3 = 1/V3- 

Saulson's 1983 numerical estimates |^ of the seismic 
gravity-gradient noise were based on seismic noise lev- 
els W{f) = 0.5 X 10~^(10Hz//)2cm/VH^ for "average 
sites" and a factor 10 lower than this for "quiet sites". 
The resulting gravity-gradient noise x{f) — T{f)W{f) 
was substantially below the projected vibrational seis- 
mic noise in (seismically well isolated) "advanced" LIGO 
interferometers [|]. 

In updating these estimates for his recent mono- 
graph |1^, Saulson noted that his original "average" 
and "quiet" sites were based on measurements at un- 
derground seismological stations. Surface sites, such 
as those chosen for LIGO and VIRGO, are far noisier 
than underground sites in the relevant frequency band, 
3Hz < / < 30 Hz, because of surface seismic waves. 
More specifically, even though the chosen LIGO sites 
(at Hanford, Washington and Livingston, Louisiana) are 
among the more quiet locations that were studied in the 
LIGO site survey, their noise at typical times is approxi- 
mately isotropic [Z{f) ~ X{f) ^ W{f)] and has approx- 
imately the following form and magnitude [|l^,0 

~ cm 

W{f) = 1x10-'^^= atl</<10Hz, 
VHz 



1 X 10 



_^ cm /lOHzV 

vTz \~r) 



at / > 10 Hz. (1.27) 



This so-called standard LIGO seismic spectrum is 20 
times larger than at Saulson's original "average" sites for 
/ > 10 Hz. Correspondingly, Saulson pointed out in his 
update, the seismic gravity-gradient noise may stick up 
above the vibrational seismic noise in "advanced" LIGO 
interferometers]^ On the other hand, at very quiet times 
— at night and with winds below 5 mph — the LIGO 
seismi c gro und noise W{f) can be as low as ~ 1/10 the 
level ( 1.27 ), thereby pushing Saulson's seismic gravity- 
gradient noise well below the vibrational seismic noise of 
an "advanced" LIGO interferometer. 



t Saulson informs us that in evaluating the noise at the LIGO 
sites, he made an error of \/3; his transfer function and the 
standard LIGO seismic spectrum actually predict a noise level 
y/S smaller than he shows in Fig. 8.7 of his book [^. When 
this is corrected, his predicted noise, like ours, is below the 
"advanced" LIGO noise curve, though only slightly so near 
10 Hz. 



TABLE I. Reduced transfer functions /3' predicted for Han- 
ford and Livingston by our 4-layer models; and /3^, the 
value of /3' when enough Love waves are added to bring the 
anisotropy ratio down to the quiet-time values observed at 
the two sites (^ ~ 1 for Hanford, A — 0.6 for Livingston). 



Modes 


Hanford 


Hanford 


Livingston 


Livingston 




P' 








RF / < lOHz 


0.4-0.85 


0.35-0.6 


0.65-0.9 


0.35-0.45 


RF / > lOHz 


0.85 


0.6 


0.65-0.9 


0.35-0.45 


RS 


0.4-1.4 


0.4-1.05 


0-1.2 


0-0.9 


RP 


0-0.15 


0-0.15 


0.02-0.13 


0.01-0.06 



D. Our analysis and transfer function 

Saulson's new, more pessimistic estimates of the seis- 
mic gravity gradient noise triggered us to revisit his 
derivation of the transfer function T(/) from seismic 
groimd motions to detector noise. Our analysis consists 
of: 

(i) splitting the ambient seismic motions into Love and 
Rayleigh modes (body of this paper and appendices); 

(ii) computing the reduced transfer function for each 
mode and for models of the geological strata at each 
LIGO site (body and appendices); 

(iii) using seismic measurements at the LIGO sites and 
geophysical lore based on other sites to estimate the mode 
mixture present at the two sites under both quiet and 
noisy conditions (this section); and 

(iv) evaluating for these mode mixtures the expected 
reduced transfer function and resulting noise (this sec- 
tion). 



1. Our reduced transfer functions 

Table ^ summarizes the results of our model computa- 
tions for each LIGO site. Shown there are the range of 
computed reduced transfer functions for specific types 
of Rayleigh modes, and the range of net reduced trans- 
fer functions that would result if each Rayleigh mode 
were mixed with enough Love waves to bring its (often 
rather high) anisotropy ratio A down to the level typical 
of quiet times at the LIGO sites (y^ ~ 1.0 at Hanford 
[|l||, ^ ~ 0.6 at Livingston jll). 

The modes shown in Table| are the RF mode, the RS 
modes with no sign of RP admixture, and the RP modes. 
The RF and RS modes usually have /?' in the range 0.4 
to 1.2, though in special cases it can sink toward zero. 
By contrast, the RP modes always have small /3': to 
0.15. This marked difference arises from the fact that for 
RF and RS the (largely S-wave) surface source tends to 
dominate over the (entirely P-wave) subsurface source; 
while for RP, mass conservation guarantees that the two 
sources (both largely P-wave) will be nearly equal, but 
opposite in sign, and will nearly cancel. (If the surface 
source were absent, the pattern would be reversed: the 



7 



subsurface source J TZ, which arises from compressional 
density perturbations, tends to be weak for RS modes 
because they consist primarily of non-compressional S- 
waves, but is strong for RP modes since they consist pri- 
marily of compressional P- waves; so /3' would be small 
for RS and large for RP.) 



2. Modes actually present and resulting seismic noise 

There is little direct evidence regarding which modes 
contribute to the ambient surface motions and thence to 
the gravity-gradient noise at the LIGO sites during quiet 
times. Past seismic measurements do not shed much light 
on this issue. In the concluding section of this paper (Sec. 
0), we shall propose measurements that could do so. 

Fortunately, the nature of the ambiently excited modes 
has been studied at other, geophysically similar sites 
(horizontally stratified alluvia over bedrock). The pre- 
ponderance of evidence suggests that at quiet times the 
surface motions at such sites and in our frequency band 
are due to a mixture of Love waves and the fundamen- 
tal Rayleigh mode RF plus perhaps a few low order RS 
modes |15Hl9|]. In at least one case, some amount of 
RP excitation is also seen these RP excitations are 
ascribed to "cultural noise" (noise generated by human 
activity of some sort) near the measurement site. Deep 
borehole measurements indicate that RP dominates at 
very great depths (^ 5 km) ; this is probably not rel- 
evant to our analysis, however. It merely indicates that 
very deep down, the majority of the surface waves have 
damped away, leaving only some residual RP modes. The 
deep motions are typically an order of ma gnit ude or two 



smaller than the surface motions; cf. Sec. IC 



On this basis, we presume that at quiet times the net 
reduced transfer function is about that for the RF mode, 
with enough admixed Love waves to bring the net A down 
to the typical quiet-time values of 1.0 for Hanford and 0.6 
for Livingston. In other words, /J^et is about equal to 
for the RF mode: 



net, quiet times 



0.35-0.45 at Livingston, 
0.35-0.6 at Hanford. (1.28) 



We have folded these quiet-time estimates for f3' into 



the standard LIGO seismic spectrum (1.27) to obtain the 



gravity-gradient noise estimates shown as the dark gray 
band in Fig. 2. The thickness of the band indicates the 
range of our /3' [Eq. ([1^)1^.35 to 0.6. To produce this 
plot, we took expression ( |l.8| ) for the transfer function 
T(/) with 7 and F set to unity, so /3 = /3' [cf. Eq. (|1.17D]. 
Then, we multi plied this by the standard LIGO seismic 
spectrum (1.27) for the ground displacement with an as- 
sumed density p ~ 1.8 g/cm^. This yields 



tSGG 



(/) 



^' 6 X 10-23 /lOHz\ 



0.6 



Hz 



3Hz < / < 10 Hz, 




10-23 _ 



7 8 9 10 20 30 

frequency, Hz 

FIG. 2. Seismic gravity-gradient noise in a LIGO interfer- 
ometer. In this figure, we assume that the direction-averaged 
spectrum of earth displacemen ts ha s the form of the standard 
LIGO seismic spectrum, Eq. (1.27). The edges of the gray 
bands are for the indicated values of the reduced transfer 
function (assumed equal to /3; i.e., for 7 and F approxi- 
mated as unity). The dark gray band is our estimate of the 
range of noise for quiet times. The gray bands, both light 
and dark, are for noisy times, assuming the standard LIGO 
seismic spectrum ( |l.2'i| ). At very quiet times, the ground 
spectrum can be a factor ~ 10 smaller than (1.27), which 



will lower these bands accordingly. Conversely, at noisy times 
the ground spectrum can be larger, raising these bands. Also 
shown for comparison is the projected noise in an "advanced" 
LIGO interferometer, and the standard quantum limit (SQL) 
for an interferometer with one tonne test masses. The SQL 
is the square root of Eq. (122) of Ref. The "advanced" 

interferometer noise is taken from Fig. 7 of Ref. Q , with cor- 
rection of a factor 3 error in the suspension thermal segment 
(Fig. 7 of Ref. is a factor 3 too small, but Fig. 10 of that 
reference is correct, for the parameters listed at the end of the 
section "LIGO Interferometers and Their Noise"). 



/3' 6 X 10-23 ^lOHzy 



0.6 



10 Hz < / < 30 Hz, 
(1.29) 



which we plotted for the indicated values of (3' . 

At very quiet times, the ambient seismic spectrum near 
10 Hz can be as much as a factor ~ 10 lower than the 



standard LIGO spectrum assumed in Eq. (1.29) and Fig. 
2, and correspondingly the quiet-time gravity gradient 
noise can be a factor ~ 10 lower. 

At noisier times, there appear to be excitations of a 
variety of RF, RS and RP modes. For example, at the 
LIGO sites, time delays in correlations between surface 
motions at the corner and the end stations reveal hori- 
zontal propagation speeds ch 5000 m/s, corresponding 
to deeply seated RP-modes (although for the most part 
these modes are seen at frequencies too low to be of in- 
terest in this analysis — / ^ 0.2 Hz p3| , p4[ ). Moreover, 
the measured anisotropy ratios can fluctuate wildly from 
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^ 0.2 to 5 at noisy times, suggesting a wildly fluctuat- 
ing mixture of RF, RS, RP, and Love modes. Scrutinizing 
not only Table | but also the range of /?' shown in Figs. 
0, ^, |l^ and which underlie that table, and keeping 
in mind that Love modes with vanishing /?' will also be 
present, we estimate that the fluctuations of /3' at noisy 
times will be confined to the range 



net, noisy tunes 



0.15-1.4. 



(L30) 



We have folded this estimate into the standard LIGO 
seismic spectrum to obtain the upper and lower edges 
of the light gray band in Fig. 2. The grey bands, light 
and grey taken together, are our best estimate of the 
range of seismic gravity-gradient noise at noisy times, 
assuming the standard LIGO seismic spectrum. Since, 
at noisy times, the seismic spectrum can be somewhat 
higher than the standard one, the gravity-gradient noise 
will be correspondingly higher. 

For the next few years, the most important applica- 
tion of these estimates is as a guide for the development 
of seismic isolation systems and suspension systems for 
LIGO. There is not much point in pushing such systems 
so hard that the vibrational seismic or the suspension 
thermal noise is driven far below our lowest estimates of 
the seismic gravity-gradient noise [bottom of the black 
line in Fig. 2, lowered by the amount that the actual 
very quiet time spectrum falls below the standard LIGO 
spectrum (1.27)] — unless corresponding steps are taken 
to mitigate the seismic gravity gradient noise; see Sec. 

In Fig. 2 we compare our predicted seismic gravity gra- 
dient noise to the projected noise in "advanced" LIGO 
interferometers and to the standard quantum limit for an 
interferometer with one tonne test masses ("SQL"). No- 
tice that our lower bound on the seismic gravity-gradient 
noise is everywhere smaller than the "advanced" inter- 
ferometer noise, but it is larger than the SQL at frequen- 
cies below ~ 20 Hz. Our lower bound rises large enough 
below ~ 10 Hz to place limits on seismic-isolation and 
suspension-noise R&D that one might contemplate do- 
ing at such frequencies. 

The remainder of this paper is organized as follows: we 
begin in Sec. |l^ by discussing Rayleigh waves and seismic 
gravity-gradient noise in the idealized case of a homoge- 
neous half space (not a bad idealization for some regions 
of some modes at Hanford and Livingston). Then we de- 
velop multilayer geophysical models for Hanford and Liv- 
ingston and use them to derive the reduced transfer func- 
tions for the various Rayleigh modes (Sees. IH and IV). 
We conclude in Sec. ^ with a discussion of the uncertain- 
ties in our analysis and research that could be undertaken 
to reduce the uncertainties, and also a discussion of the 
physical interaction of the seismic waves with the foun- 
dations of the LIGO facilities, and of ways to somewhat 
reduce the gravity gradient noise if it ever becomes a 
serious problem in LIGO interferometers. Mathematical 
details of our analysis are confined to Appendices. Those 



Appendices may form a useful foundation for analyses of 
seismic gravity-gradient noise at other sites. 



II. HOMOGENEOUS HALF SPACE 

A. Fundamental Rayleigh mode 

As a first rough guide to seismic gravity-gradient noise, 
we idealize the LIGO sites as a homogeneous half space 
with density p, Poisson ratio v, S-wave speed cs and P- 
wave speed cp given by 

p=1.8g/cm^, ^ = 0.33, CP = 440 m/s, cs = 220 m/s. 

(2.1) 

(These are the measured parameters of the surface mate- 
rial at Livingston; for H anford , the parameters are only 
a little different; cf. Sec. [II A below.) 

This homogeneous half space can only support the RF 
mode, as mentioned in the Introduction. The theory of 
the RF mode and the seismic gravity^radient noise that 
it produces is sketched in Appendix ^ Here we summa- 
rize the results. 

The RF mode propagates with a horizontal speed ch 
that depends solely on the Poisson ratio. It is a bit slower 
than the speed of S- waves, and is much slower than P- 
waves. For the above parameters. 



CH = 0.93cs = 205 m/s ; 



(2.2) 



cf. Eq. (Bl). Correspondingly, the waves' horizontal 
wave number k and horizontal reduced wavelength are 



A 1 /lOHz 

t:— = T = -J-om — n — 

2tt k \ f 



(2.3) 



Because ch < cs < cp, RF waves are evanescent ver- 
tically: the P-waves die out with depth z as e""^*^^, and 



the SV-waves as 



where 



g = Vl - [ch/cpY = 0.88 , 
s = Vl - {ch/csY = 0.36 . 



(2.4) 



Thus, the vertical e-folding lengths for compression 
(which produces seismic gravity gradients) and shear 
(which does not) are 



Zp = — = 3.7 m 
qk 

Zs^\ = 9.0m 
sk 



10 Hz 
/lOHz 



(2.5) 



These RF waves produce substantially larger vertical 
motions than horizontal at the earth's surface. For waves 
that are horizontally isotropic, the anisotropy ratio is 



-^=^ l-f.^-2g. ='-' 



(2.6) 
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This large ratio is indicative of the fact that RF waves 
contain a large component of P-waves. As mentioned 
in the Introduction, this is substantially larger than the 
values typically observed at the LIGO sites in the band 
3Hz < / <_^30Hz — seismic measurements taken at 
those sites [ |l3yi4| show that, at quiet times, ^ ~ 1.0 at 
Hanford, A ~ 0.6 at Livingston. Thus, RF waves can- 
not alone be responsible for the seismic motions. To the 
extent that our homogeneous-half-space model is realis- 
tic, RF waves must be augmented by a large amount of 
horizontally-polarized S- waves ("SH- waves"), which have 
A = 0. 

RF waves produce a reduced transfer function 



3(1 + s2 - 2qf 



2(l + 52)[(l + .s2)(l+q2)_4^s] 



0.86. (2.7) 



This /?' is produced primarily by the surface source in 
Eq. Ql 24| ); if there were no surface source, the subsurface 
term J TZ (arising solely from the P-wave compressions) 
would produce the far smaller value f3' = 0.17. When 
the RF waves are augmented by enough Love waves to 
reduce the net to 1.0 (Hanford) or 0.6 (Livingston), 
they produce a net reduced transfer function [Eq. ( l.K: ) 
with primes added to the /3's] 



= 0.59 (Hanford), 0.40 (Livingston). 



(2.8) 



As we shall see in the next two sections, the earth is 
strongly stratified over the relevant vertical length scales 
at both Hanford and Livingston, and this gives rise to 
significant differences from the homogeneous-half-space 
model. Nevertheless, as discussed in the Introduction 



(Sec. ID 2 ), it is likely that at quiet times the RF mode 
produces the dominant gravity-gradient noise. Stratifica- 
tion modifies this RF mode somewhat from the descrip- 
tion given here; however, as we shall see (Figs. |^ and 
|lO|) , these modifications typically alter its anisotropy ra- 
tio and reduced transfer function by only a few tens of 
percent. Thus, the homogeneous-half-space model may 
be a reasonable indicator of seismic gravity-gradient noise 
in LIGO at quiet times. 



B. P-up and SV-up waves 

The principal effect of stratification is to produce a rich 
variety of normal- mode oscillations, in which mixtures of 
SV- and P-waves resonate in leaky cavities formed by the 
strata. These oscillations are Rayleigh-mode overtones, 
whose (rather complex) theory is sketched in Appendices 
^ and ^ and discussed in Sees. [H and In this subsec- 
tion we will momentarily ignore that fact, and will seek 
insight from a much simpler analysis that gives gives re- 
sults which agree approximately, and in some cases quite 
well, with those of the Rayleigh-overtone theory. 

If the top layer (labeled by a subscript 1) has a thick- 
ness Di larger than half a vertical wavelength of the 



\ 1 

\ ; 

p\ ' 

\ / 

\^ / 
\. / 


\ 

\ 

\ 

\ 

\ 

\ 

\ 


/ 

/sv 





0.2 



0.8 



1 



0.4 0.6 

cP/cH or cS/cH 

FIG. 3. Anisotropy ratio for waves that propagate upward 
in a homogeneous half space, reflect off the Earth's surface, 
and propagate back downward. The curve "P" is for the case 
when the upward propagating waves are pure P (P-up waves), 
in which case the abscissa is cp/cn = sin ap; "SV" is for 
SV-up waves, with abscissa cs/ch = sin as . It is assumed 
that cp = 2cs; this is approximately the case for the surface 
layers at Hanford and Livingston. 



waves' oscillations, Di > (cpi/2/)/^/l — (cpi/chY for 
P-waves and similarly for S- waves [cf. Eq. (3^) below 
and associated discussion], then the trapped modes can 
be thought of as propagating upward through the top 
layer, reflecting at the earth's surface, and then prop- 
agating back downward. By ignoring the effects of the 
interfaces below, these waves can be idealized as traveling 
in a homogeneous half space. 

The behavior of these waves depends on the mixture of 
P- and SV-waves that composes them as they propagate 
upward. Because these two components will superpose 
linearly, we can decompose the mixture and treat the P- 
wave parts and SV-wave parts separately. We will call 
these components P-up and SV-up waves. In Appendix 

we derive simple analytic formulae for the anisotropy 
ratio A and reduced transfer function /3' for P-up and SV- 
up waves, and in Figs. ^ and ^ we graph those formulae. 
In these plots, for concreteness, we have chosen cs = 
cp/2. 

Consider, first, the P-up waves (so lid c urves in Figs. 
I and|). Due to Snell's law [cf. Eq. below], these 

waves propagate at an angle ap — arcsin(cp/c/f ) to the 
vertical. Such propagating waves can therefore exist only 
for ch > Cp] when ch < cp, P-waves are evanescent. For 
this reason, in the Figures we plot on the abscissa the ra- 
tio cp/cH running from to 1. When P-up waves hit the 
surface, some of their energy is converted into SV-waves 
propagating downward at an angle as = arcsin(cs/c//); 
the rest of the energy goes into reflected P-waves. The 
resulting combination of upgoing P- and downgoing P- 
and SV-waves gives rise to the anisotropy and reduced 
transfer functions shown in the figures. 

For Ch S> Cp the waves travel nearly vertically. Their 
P-components produce vertical motions, while the much 
weaker SV-waves created on refiection produce horizon- 
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0.2 0.4 0.6 0.8 1 

cP/cH or cS/cH 

FIG. 4. Reduced transfer function for P-up and SV-up 
waves in a homogeneous half space with cp — 2cg. Nota- 
tion is as in Fig. H. 



tal motions. As a result, A is large, diverging in the 
limit c/f — > oo, and decreasing gradually to near unity as 
ch Cp- As we shall see below, this is typical: when 
P-waves predominate in a wave mixture (RP modes), A 
is typically somewhat larger than unity. 

For these P-up waves, the gravity gradients produced 
by the surface source cancel those from the subsurface 
source in the limit ch 3> cp, causing /?' to vanish. As ch 
is reduced (moving rightward in Fig. the cancellation 
becomes imperfect and P' grows, though never to as large 
a value as /?' would have in the absence of the surface 
term (~ 1.3-2.4). The surface-subsurface cancellation is 
easily understood. In the limit ch ^ cp, the P-waves 
propagate nearly vertically, with vertical reduced wave- 
length for their density oscillations, 1/fcy = cp/w, that is 
small compared to the gravity-gradient e-folding length 
-2^sgg = 1/k = ch/oj, over which the waves' sources 



are integrated in Eq. (1.24) to produce the gravitational 



force. Therefore, the gravity gradients come from many 
vertical wavelengths, with adjacent ones weighted nearly 
equally. Because of mass conservation, the surface source 
plus the top quarter wavelength of subsurface source 
(multiplied by p) constitute the mass per unit area that 
has been raised above a node of the mode's displacement 
eigenfunction; and correspondingly their sum vanishes. 
Below that node, alternate half wavelengths of the sub- 
surface source cancel each other in a manner that gets 
weighted exponentially with depth, e"*^^; their cancella- 
tion is excellent in the limit 1/fcy ^ 1/A;, ie., ch S> cp. 

Turn now to the SV-up waves. Upon reflection from 
the surface, these produce a mixture of downgoing SV- 
and P-waves. This mixture gives rise to the anisotropy 
and reduced transfer functions shown dashed in Figs, y 
and|^. Again by Snell's law, SV-up waves propagate at an 
angle as ~ &rcs\n{cs / ch) to the vertical; thus, propaga- 
tion is possible only for ch > cs, and so we plot on the ab- 
scissa cs/ch running from to 1. When ch > '^cs — cp 
(left half of graphs), the downgoing P-waves generated at 
the surface can propagate; when ch < '2cs = cp (right 
half of graphs), the downgoing P-waves have imaginary 



propagation angle a p and thus are evanescent (decay ex- 
ponentially with depth). This is analogous to the phe- 
nomenon of total internal reflection which one encoun- 
ters in elementary optics. The downgoing P-waves are 
the sole subsurface source of gravity-gradient noise, and 
since they are only a modest component of the SV-Up 
mode, the subsurface source is small. The SV-waves pro- 
duce no subsurface source (no compressions), but they 
produce a large surface source (large surface vertical mo- 
tions) . This surface source is the dominant cause of the 
gravity-gradient noise and predominantly responsible for 
the rather large reduced transfer function shown in Fig. 
^. Note that the maximum value, (3' ~ 1.4, is the same 
as the largest /?' for RS modes in our 4-layer models of 
the LIGO sites (Table |). 

When propagating more or less vertically (ch > 2C5), 
these SV-up waves produce small anisotropics {A < 0.4 
— large horizontal motions and small vertical motions). 
When they propagate more or less horizontally, A is 
large. The divergence of A at cs/ch — l/\/2 — 0.707 
(as = 7r/4) occurs because the SV-up waves at this an- 
gle generate no P-waves upon reflection; they only gen- 
erate downgoing SV-waves, and the combination of the 
equal-amplitude up and down SV-waves produces purely 
vertical motions at the earth's surface. At frequencies 
/ ^ 20 Hz, mode RSI at Hanford can be approximated 
as an SV-up mode and exhibits this behavior; cf. Sec. 
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III. HANFORD 
A. Hanford geophysical structure 

At the LIGO site near Hanford, Washington, the top 
220 m consists of a variety of alluvial layers (fluvial and 
glacio-fluvial deposits of the Pliocene, Pleistocene, and 
Holocene eras; coarse sands and gravels, fine sands, silts, 
and clays, in a variety of orders). The upper 40 m are 
dry; below about 40 m the alluvium is water-saturated. 
From the base of the alluvium (220 m) to a depth of 4 
km lies a sequence of Columbia River basalts, and below 
that, bedrock [||,|6|. 

The density of the alluvial material is p ~ 1.8 g/cm^, 
independent of layer. Velocity profiles (cp and cs as 
functions of depth z) have been measured at the site by 
contractors in connection with two projects: LIGO |35| ] 
and the Skagit nuclear power plant (which was never 
constructed). We have relied primarily on the Skagit 
report because it contains more detailed information over 
the range of depths of concern to us, and because there 
is a serious discrepancy between the two reports in the 
depth range 5-25 m, which contributes significantly to 
the seismic gravity gradients. The Skagit velocities there 
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are moreplausible than the LIGO one^. 

Table || shows velocity profiles as extracted from the 
Skagit report. Notice the overall gradual increase in both 
wave speeds. This is due to compression of the alluvia 
by the weight of overlying material, with a consequent 
increase in the areas of the contact surfaces between ad- 
jacent particles (silt, sand, or gravel) |2^. Notice also 
the sudden increase of cp and j/ at 40 m depth, due to a 
transition from dry alluvia to water-saturation; the wa- 
ter contributes to the bulk modulus but not the shear 
modulus, and thence to cp but not cg. Notice, finally, 
the large jump in both cp and cs at the 220 m deep 
transition from alluvial deposits to basalt. 

We have been warned by geophysicist and seismic en- 
gineer colleagues that we should not place great faith 
in all the details of measured velocity profiles such as 
this one; and the discrepancies between the Skagit and 
LIGO velocity-profile measurements have reinforced this 
caution. As a result, from computations based on these 
velocity profiles (and similar profiles at Livingston), we 
can only expect to learn (i) the general nature of the 
modes to be expected at each LIGO site, (ii) how those 
modes' characteristics are infiuenced by the velocity pro- 
files, (iii) the range of anisotropy ratios A and reduced 
transfer functions /?' to be expected at each site, and 
(iv) how A and depend on the velocity profiles and 
the modes' characteristics. We cannot expect the com- 
puted, mode-by-mode details of A{f) and /?'(/) to be 
accurate — except, perhaps, for the shallowly seated RF 
mode. Nevertheless, the insights that we do gain from 
such computations should be of considerable help in fu- 
ture studies of seismic gravity-gradient noise and future 
attempts (if any) to mitigate it. 

In this spirit, we have simplified our calculations by ap- 
proximating the measured Hanford velocity profiles (Ta- 
ble H) with their twelve distinct layers by the simpler 
four-layer model shown in Table [II. Layers 1 and 2 are 
dry alluvia, layer 3 is water-saturated alluvium, and layer 
4 is basalt. 



TABLE II. Velocity profiles at the Hanford LIGO site, as 
extracted from Table 2.5-3, Fig. 2.5-10, and Sec. 2.5.2.5 
of the Skagit Report |Q. These velocities are based on (i) 
cross-hole measurements (waves excited in one borehole and 
measured in another) down to 60 m depth; (ii) downhole mea- 
surements (waves excited at surface and arrivals measured in 
boreholes) from z — 60 mto z = 175 m; (iii) extrapolations of 
downhole measurements at other nearby locations, and sur- 
face refraction measurements (waves excited at surface and 
measured at surface) at the LIGO site, from z = 175 m down 
into the basalt at z > 220 m. The downhole measurements at 
one well (Rattlesnake Hills No. 1) have gone into the basalt 
to a depth of 3230 m. Depths are in meters, velocities are in 
m/s. 



Depths 


Cp 


cs 


u 


0-12 


520 


270 


0.32 


12-24 


820 


460 


0.27 


24-32 


1000 


520 


0.31 


32-40 


1260 


530 


0.39 


40-50 


1980 


560 


0.46 


50-80 


2700 


760 


0.46 


80-110 


2700 


910 


0.44 


110-160 


1800 


610 


0.44 


160-210 


2400 


910 


0.42 


210-220 


2900 


1200 


0.40 


220-250 


4900 


2700 


0.28 


250-3230 


5000-5700 competent 

basalt flows 
4000-5500 interbeds 







B. Hanford model results 

The horizontally stratified geologies at Hanford and 
Livingston support a variety of Love and Rayleigh modes. 



TABLE III. Four-layer model for the velocity proflles at 
the Hanford LIGO site. Notation: n — layer number, Z)„ 
— layer thickness, cp„ — P-wave speed in this layer, cs-n. — 
S-wave speed in this layer, Vn — Poisson ratio in this layer. 
Depths and thicknesses are in meters, speeds are in m/s. 



■''The report prepared for LIGO |3^ claims cp = 1400 m/s, 
Cs ~ 370 m/s, corresponding to a Poisson ratio of = 0.46. 
This could be appropriate for water-saturated materials at 
this depth, but is not appropriate for the dry materials that 
actually lie there. The Skagit report jsij shows two layers in 
this range of depths: one with cp = 520 m/s, cs — 270 m/s, 
for which v = 0.32; the other with cp = 820 m/s, cs = 460 
m/s, for which v = 0.27. For dry alluvia, these values are 
much more reasonable than i/ — 0.46. We thank Alan Rohay 
for bring this point to our attention. 
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Depths 


D„ 


CPn 


CSn 




1 


0-12 


12 


520 


270 


0.32 


2 


12-40 


28 


900 


500 


0.28 


3 


40-220 


180 


2400 


700 


0.45 


4 


220-4000 


3780 


4900 


2700 


0.28 
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(For the general character of Love and Rayleigh modes 
see, e.g., Refs. (H^^J] and the brief discussion in the in- 
troduction of this paper.) We shall focus on Rayleigh 
modes in this section, since they are the sole producers 
of seismic gravity-gradient noise. 

In each geological layer, consider a specific Rayleigh 
mode. It consists of a superposition of plane-fronted P- 
and SV- waves. Because each layer is idealized as ho- 
mogeneous, the mode's SV- and P-waves are decoupled 
within the layer. However, they are coupled at layer in- 
terfaces and the earth's surface by the requirement that 
material displacement and normal stress be continuous 
across the interface (or with the atmosphere in the case 
of the earth's surface) . The details of this coupling and 
its consequences are worked out in Appendix |y. 

In each layer, the mode's P- and SV-components prop- 
agate at different angles to the vertical: ap„ for the P- 
waves in layer n and asn for the SV- waves. However, 
the components must all move with the same horizontal 
speed 



CPn 



CSr, 



k 



sm ap„ smasri 



(3.1) 



(Snell's law), and they must all have the same horizontal 
wave number k and frequency uj — 2-Kf. 

Each mode can be characterized by its dispersion re- 
lation for horizontal motion Lo{k), or equivalently cnif)- 
It will be helpful, in sorting out the properties of the 
modes, to understand first what their dispersion rela- 
tions cuif) would be if their SV-wave components were 
decoupled from their P-wave components. We shall do 
so in the next subsection, and then examine the effects of 
coupling in the following subsection. Note that we shall 
ignore the effects of damping in these two subsections, 
since the lengthscales involved are less than (or at most 
of the same order as) the dissipation lengthscales of both 
P- and SV-waves [c/. Eq. (IT" 



1. P-SV decoupling approximation 

Recall that we denote by RPn the nth Rayleigh mode 
of P-type and by RSn the nth Rayleigh mode of SV-type. 
In the approximation of P-SV decoupling, Mode RPn 
with horizontal speed ch propagates from the earth's 
surface through sequences of strata (generating no SV- 
waves) until it reaches a depth Dp where cp first exceeds 
Ch- At that location, it reflects and returns to the sur- 
face, and then is reflected back downward. The mode's 
dispersion relation cnif) is determined by the resonance 
condition that the reflected waves arrive at the surface 
in phase with the original downgoing waves. 

This resonance condition is evaluated most easily by 
following the (locally) planar waves vertically downward 
to their reflection point (the location z = Dp where cp 
first reaches ch) and then back up, thereby returning 
precisely to the starting point. On this path, the vertical 
component of the wave vector is 
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CO 650 
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FIG. 5. Dispersion relations for the 4-layer Hanford 
model, as computed using the P-SV decoupling approxima- 
tion, Eqs. (O) and (E|). 



fcy = — cosap = — a/1 — {cp/chY = kcotap , (3.2) 

Cp Cp 



where we have used Snell's law ( |3.l| ) to infer cos ap = 
y/l — [cp/chY- The waves' corresponding waves' total 
round-trip phase shift is 



LO 



A$ = 2 / —^l^{cp/cHfdz + (5$i„tcrfacc« ■ (3.3) 
Cp 

Here (5$interfaces IS the total phase shift acquired at the in- 
terfaces between strata and upon reflecting at the earth's 
surface. Setting oj — 27r/ and imposing the resonance 
condition A<I> — 2m:, we obtain the following dispersion 
relation for mode RPn: 



/ 



n - ((5$intcrfaccs/27r) 



(3.4) 



Similarly, for mode RSn the dispersion relation is 

n - ((5<I'intcrfaccs/27r) 



Ds 



(3.5) 



where Ds is the depth at which cs first reaches ch. 

Figure [s] shows these decoupling-approximation disper- 
sion relations for our 4-layer model of cp{z) and cs{z) 
(Table III). For the RS- waves, the total interface phase 
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shift has been set to i^^'intcrfaccs = i", which would be the 
value for a single layer with a huge rise of cs at its base. 
For the sole RP mode shown, RPl, it has been set to 
'^'I'intorfaccs — 7''/2, which is a fit to the dispersion rela- 
tion with P-SV coupling (Fig. |[ to be discussed below). 

Notice that for fixed horizontal speed ch, the lowest 
RP mode, RPl, occurs at a much higher frequency / 
than the lowest RS mode, RSI. This is because of the 
disparity in propagation speeds, cp = severalxcs. Notice 
also the long, flat plateaus in cnif) near ch = cs2 = 500 
m/s and especially C53 — 700 m/s for the RSrt modes, 
and near ch = cp2 = 900 m/s and ch = cp3 = 2400 m/s 
for RPl. Mathematically these are caused by the van- 
ishing squ are roots i n th e denominators of the dispersion 
relations (3^) and (3^). Physically they arise because 
the mode's waves "like" to propagate horizontally in their 
deepest layer. At high frequencies {e.g., / > 10 Hz for 
CH — CS3 = 700 m/s), several modes propagate together 
nearly horizontally in that deepest layer. 



2. Effects of P-SV coupling on dispersion relations 

Figure ^ shows the dispersion relations cuif ) for the 
lowest 8 modes of our 4-layer model at Hanford, with 
P-SV coupling included. These dispersion relations were 
computed using the multilayer equations of Appendix |^. 
We shall now discuss these various dispersion relations, 
beginning with that for the fundamental mode, which is 
labeled RF in the figure. 

Mode RF was studied in Sec. 
mogeneous half space. It is vertically evanescent in both 
its P- and SV-components (except at low frequencies in 



II A for an idealized ho- 



the top layer); for this reason, it did not show up in 
our idealized decoupling-approximation dispersion rela- 
tion (Fig. ^). At frequencies / ^ 10 Hz, its vertical 
e-folding lengths Zp and Zs [Eqs. ( p.5D ] are both short 
enough that it hardly feels the interface between layers 
1 and 2, and the homogeneous-half-space description is 
rather good. Below 10 Hz, interaction with the interface 
and with layer 2 pushes ch up. 

By contrast with the P-SV-decoupled Fig^^ ||, every 
Rayleigh overtone mode RPn or RSn in Fig. |6| now con- 
tains a mixture of SV- and P-waves. This mixture varies 
with depth in the strata and is generated by the same 
kind of i nterface reflection and refraction as we met in 
Sec. [IB for SV-up and P-up waves. In most regions of 
the (c/f,/) plane, the mode mixtures are dominated ei- 
ther by SV- or P-waves — the ratio of energy in one wave 
type to that in the other is > 2. 

In the vicinity of the wide gray band marked RPl, the 
modes are predominately of RP type; away from that 
vicinity they are predominately RS. The location of the 
RPl band has been inferred from the computed S- and 
P-wave amplitudes. Notice how well it agrees with the 
decoupling approximation's RPl dispersion relation (Fig. 
Away from the RPl band, the dispersion relation 
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FIG. 6. Dispersion relations for the 4-layer Hanford model, 
including coupling between P- and SV-waves produced at 
boundaries between layers and at the earth's surface. 



for each RSn mode is reasonably close to its decoupling- 
approximation form (compare Figs. ^ and |^). As each 
mode nears and crosses the RPl band, its dispersion re- 
lation is distorted to approximately coincide, for awhile, 
with the RPl shape. Correspondingly, all its other prop- 
erties become, for awhile, those of an RP mode. 



3. Anisotropy ratios and reduced transfer functions 

Figure ^ shows the anisotropy ratio A and reduced 
transfer function (3' for the lowest eight modes of our 
4-layer model of Hanford. These were computed using 
the multilayer equations of Appendix |c[ with dissipation 
neglected. On the Figure, the mode names "RSn" have 
been shortened to "n" , and "RF" to "F" . The bottom set 
of graphs is the value that the net reduced transfer 
function would have if the mode of interest were mixed 
with enough Love waves to reduce the net anisotropy 
ratio to the value -4nct — 1.0 typical of measured seismic 
spectra at Hanford during quiet times . 

Fundamental Mode RF. Above 10 Hz, mode RF has 
A ~ 2.2, 13' ~ 0.84, and (3'^ ~ 0.58, in accord with our 
homogeneous-half-space model (Sec. ^. Below 10 Hz, 
coupling of the RF mode to layer 2 produces a growth 
of the subsurface source to partially cancel the surface 
source, and a resulting fall of /?' to 0.4 and to 0.35. 
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FIG. 7. Properties of the lowest 8 modes of the 4-layer 
Hanford model, including coupling between P- and SV-waves 
produced at boundaries between layers and at the earth's sur- 
face. 



RS Overtones. In RS regions (away from the RPl 
band) the overtone modes RSn generally have A ^ 1 so 
~ (3' — little or no admixed Love waves are needed to 
bring the anisotropy down to 1.0. The value of (3' ranges 
from ^ 0.4 to 1.4 in the RS regions; but when the RPl 
mode is nearby in the CH~f plane, its admixture drives 
/?' down to < 0.2. 

Mode RSI shows characteristic "SV-up" behavior near 
25 Hz (compare Fig. ^ with Figs. ^ and ^). Its A has a 
very large resonance and its subsurface source (not shown 
in the figures) has a sharp dip to nearly zero, resulting 
from 45° upward propagation of its SV-component in the 
top layer and no production of P-waves upon reflection. 
At frequencies above our range of interest, this same SV- 
Up behavior will occur in successively higher RSn modes. 

RPl Mode. The region of RPl behavior is shown as 
thick gray bands in Fig. (c/. the bands in Fig. The 
RPl reduced transfer function is small, < 0.15, due to 
the same near-cancellation of its surface and subsurface 



sources as we met for P-Up waves in Sec. II B and Fig. IJ. 
As each RS mode crosses the core of the RPl region, its 
(3' shows a dip and its anisotropy shows a peak, revealing 
the temporary transition to RP behavior. 

Higher-order RP Modes. The higher-order RP 
modes {n = 2, 3, ...) in our frequency band will lie in the 
vicinity of RSn overtones with n > 8. We expect these 
RPn modes to show similarly small reduced transfer func- 
tions to those for RPl, but we have not attempted to 
compute them, with one important exception: high-order 
RP modes that travel nearly horizontally in Hanford's 
^ 4 km thick basalt layer. We consider these modes in 
the next subsection. 



4- RP modes that travel horizontally in the basalt 



As discussed in the Introduction (Sec. ID 2), the 
ground motions at the Hanford corner and end stations 
sometimes show time delays in correlated motion, cor- 
responding to wave propagation speeds of ^ 5000-6000 
m/s |l4|j2^. These motions must be due to wave modes 
that travel nearly horizontally in the ~ 4 km thick basalt 
layer at the base of the alluvium, or in the bedrock be- 
neath the basalt. We have computed the properties of 
such wave modes for the case of horizontal propagation in 
the basalt layer — layer 4 of our 4-layer Hanford model. 

Because of the many closely spaced modes in the 
relevant (ch,/) region {ch a little larger than cp4 = 
4900m/s, 3Hz < / < 30 Hz), it is not reasonable, or even 
of interest, to compute their dispersion relations explic- 
itly. Instead, we have assumed an idealized dispersion 
relation ch — 4910 m/s independent of frequency. 

The basalt layer is so thick that nearly horizontally 
propagating waves will be substantially damped in trav- 
eling from its lower face to its upper face and back; and, 
the S-waves will be much more strongly damped than 
the P-waves. For this reason, we idealize these waves as 
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FIG. 8. Properties of RP modes that propagate nearly 
horizontally in layer 4 (~ 4 km thick basalt layer) at Hanford 
{ch slightly larger than cp4 = 4900 m/s), including the effects 
of dissipation in the alluvium above the basalt. 



purely P-up as they impinge from the basalt layer 4 onto 
the layer 3-4 interface. These P-up waves at interface 
3-4 are treated as a source for other wave components in 
all 4 layers. 



For these waves, dissipation [Eqs. (O) and associated 
discussion] may be more important than for the RF, RS 
and RPl modes, which were treated above as dissipation- 
less. We therefore include it in our analysis. We do so in 
the 4-layer equations of Appendix ^ by giving the sound 
speeds appropriate imaginary parts. 



3?(cp„) ~ 2Qp 

S(cg„) 1_ 

3?(cs„) ~ 2Qs 



= -0.015 , 
= -0.03 , 



(3.6) 



while keeping their real parts equal to the values shown in 
Table III . We have solved the resulting multilayer equa- 



tions numerically, obtaining the anisotropy ratios and re- 
duced transfer functions shown in Fig. ^. 

The peaks in at / ~ 3, 11, 19, and 27 Hz [frequency 
separation A/ = cp^/ ^21)3 -\/l — Cpg/c^^ = 8 Hz] are 

associated with resonant P-wave excitations of layer 3 
and their influence on layer 2 and thence on layer 1; cf . 
the decoupling-approximation dispersion relation (B.4). 
The slightly smaller peaks at / ~ 6, 14, and 23 Hz are due 
to resonant S-wave excitations of layer 2. The oscillations 
in both A and (less obviously) f3' with frequency A/ = 

CP2/ ^21)21/1 — Cp2/c|f) = 16 Hz are associated with 
resonant P-wave excitations of layer 2. 

Despite the complexity of these deeply seated RP 
modes, with various types of resonant excitations of var- 
ious layers, and despite the fact that the seismic gravity 
gradients arise from depths Zggg = 1/A: '--^ 25 to 250 m 
so great that the top three layers all contribute, these 
modes exhibit the same range of values of /3' as mode 
RPl: to 0.15. Here, as for RPl and for P-up waves 
in a homogeneous half space, /?' is small because of near 
cancellation of the gravity gradients from the P-wave sur- 



face and subsurface sources. 



5. Summary of Hanford model results 

The most important of the above results are those for 
the reduced transfer functions of the various modes at 
Hanford. They are summarized in Table ^ and their im- 
plications are discussed in the Introduction, Sec. [D2. 



IV. LIVINGSTON 
A. Livingston geophysical structure 

At the LIGO site near Livingston, Louisiana, the ge- 
ological strata consist of alluvial deposits laid down by 
water flowing into the Gulf of Mexico. As the ocean level 
has risen and fallen, alluvial terraces of varying thickness 
have been formed. This alluvium (layers of clay, silt, 
sand, and gravel in various orders) is of the Holocene, 
Pleistocene, and Phocene eras going down to a depth 
of about 700 m, and compacted alluvium of the upper 
Miocene and earlier eras below that. These sedimentary 
deposits extend down to a depth of about 3 km ||2^] be- 
fore reaching bedrock. 

For our analysis the principal issue is the vertical ve- 
locity profiles cp{z) and cs{z). The primary difference 
between Livingston and Hanford is the depth of the wa- 
ter table: it is only about 2 m down at Livingston, versus 
about 40 m at Hanford. This difference should cause cp 
to soar to about 1600 m/s at depths of a few meters at 
Livingston; it only does so roughly 40 m down at Han- 
ford. 

The only measurements of the Livingston velocity pro- 
files that we have been able to find are those performed 
in a site survey for LIGO . Those measurements only 
include eg, not cp, and only go down to a depth of 15 m. 
Accordingly, we have had to estimate the velocity profiles 
from these sparse data and from the lore accumulated by 
the geophysics and seismic engineering communities. 

That lore suggests that cs should increase as about 
the 1/4 power of depth |26j . (This increase is due to the 
fact that the shear restoring force must be carried by the 
small-area interfaces between the grains of gravel, sand, 
silt, or clay; the weight of overlying material compacts 
the grains, increasing the areas of their interfaces.) We 
have fit the measured cs{z) in the top 15 m ({7 ft, 700 
ft/s}, {21 ft, 810 ft/s}, {50 ft, 960 ft/s}) (Ref. Ap- 
pendix B, plate 7) to a 1/4 power law, adjusting the fit 
somewhat to give speeds at greater depths in rough ac- 
cord with measurements at a similar sedimentary site in 
Tennessee S]. Our resulting fit is 



cs = 185m/s(l + z/2.9m)i/'' 



(4.1a) 



A combination of theory and phenomenology [Eqs. 
(6.24), (6.26) of [p6l and associated discussion] tells us 
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TABLE IV. Four-layer model for the velocity profiles at 
the Livingston LIGO site. Notation and units are as in Table 

inl 



1000 



n 


Depths 




CPn 


CSn 




1 


0-5 


5 


440 


220 


0.33 


2 


5-105 


100 


1660 


400 


0.47 


3 


105-905 


800 


1700 


700 


0.40 


4 


905-3005 


2100 


1900 


1000 


0.31 



that in these water-saturated aUuvia, the material's Pois- 
son ratio should be about 



0.39 f^^V 



(4.1b) 



(The Poisson ratio goes down gradually with increasing 
compaction and increasing cs because water is playing a 
decreasing role compared to the grains.) The standard 
relation 



Cp = cs 



l-2v 



(4.1c) 



combined with Eqs. (4.1a) and (4.1b), then gives us the 
vertical profile for cp . 

These profiles are valid only in the water-saturated re- 
gion. Although the water table is at ~ 2 m, measure- 
ments elsewhere suggest that one may have to go 
downward an additional several meters before the effects 
of the water on cp will be fully felt. Accordingly, we ex- 
pect Cp ~ 2cs in the top ~ 5 m at Livingston, fo llowe d 
by a sharp rise to the values dictated by Eqs. (tl.la )- 



( 4.1c ), though in our final conclusions (Sec. 1 V C 2 ) , we 
shall allow for the possibility that the sharp rise occurs 
at anywhere from 2 to 5 m depth. 



B. Livingston 4-layer model 

We have fit a four-layer model to these estimated Liv- 
ingston velocity profiles. Our fit is shown in Table 
This model is the primary foundation for our exploration 
of seismic gravity gradients at Livingston. As discussed 
above, it principally differs from the 4-layer Hanford 
model by the rapid increase of cp at 5 m depth at Liv- 
ingston, due to the higher water table. All other differ- 
ences have a much more minor influence on the seismic 
gravity-gradient noise. 



C. Livingston model results 

1. Mode overview 

Because the top, unsaturated layer is so thin, RP 
modes cannot resonate in it in our frequency band; and 
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FIG. 9. Dispersion relations for 4-layer Livingston model, 
including coupling between P- and SV-waves produced at 
boundaries between layers, for the fundamental mode and the 
lowest 10 RS modes. 



because water makes cp so large just below the top layer, 
the RP modes in our band can only propagate at a cor- 
respondingly high speed, ch > 1660 m/s. The lowest 10 
RS modes, by contrast, are confined to speeds ch ^ 1000 
m/s. As a result — in contrast to Hanford — there is 
no mixing between these lowest RS modes and the RP 
modes. The RS modes have purely RS character, with 
no significant RP admixture. 

In the next section we shall study the lowest 10 RS 
modes along with the fundamental mode. In the follow- 
ing section, we shall examine the lowest RP mode. 



2. RF and RS modes 

We have computed the dispersion relations, anisotropy 
ratios, and reduced transfer functions for modes RF and 
RSl-10 in our 4-layer Livingston model, using the multi- 
layer equations of Appendix The dispersion relations 
are shown in Fig. ||. Because of the separation in the 
(c//, /) plane of these modes from the RP mode, we ex- 
pect the P-SV decoupling approximation to work quite 
well here. Indeed, the RS modes have the form one would 
expect from the decoupling approximation [Eqs. (3.5)]. 



The anisotropy ratios and reduced transfer functions are 
shown in Fig. p^ . 

RF Mode. Because the top layer is 2.5 times thinner 
in our Livingston model than at Hanford, the frequency 
at which the RF mode becomes like that of a homoge- 
neous half space is 2.5 times higher: ~ 25 Hz compared to 
^ 10 Hz. Only above ^ 25 Hz do the mode's properties 
asymptote toward their homogencous-half-space values 
of CH = 205 m/s, A = 2.2, 13' = 0.86 and = 0.40. At 
lower frequencies, interaction with layer 2 pushes /?' into 
the range 0.65-0.9, and /9J^ into the range 0.35-0.45. 

It is possible that the effects of water saturation will 
cause Cp to shoot up at depths shallower than the 5 m 
assumed in our model; a transition anywhere in the range 
2 m ^ z 5; 5 m must be considered reasonable. If the 
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FIG. 10. Properties of the lowest 10 RS modes and the 
RF mode of the 4-layer Livingston model. (Modes RS8 and 
RS9 are not shown; their curves are sandwiched between 7 
and 10.) 



transition in fact occurs at depths shallower than 5 m, 
the peaks of /?' and will be pushed to correspondingly 
higher frequencies. Thus, we must be prepared for the 
RF mode to have (3' anywhere in the range 0.65-0.9, and 
in the range 0.35-0.45 at just about any frequency in 
our band of interest. 

RS Modes. In our frequency band, the RS modes 
have negligible excitation in layers 3 and 4, and their P- 
waves are evanescent in layers 2, 3 and 4. As a result, 
these modes can be well approximated by SV-up waves 
in layer 2, impinging on the layer 1-2 interface. We have 
verified this by computing their anisotropics and reduced 
transfer functions in this 2-layer SV-up approximation 
by the method outlined at the beginning of Appendix 
The results for A and /?' which relied on the 4-layer 
dispersion relations of Fig. 0, agree to within a few per 
cent with those of our 4-layer model (Fig. |l^) except 
at frequencies below 5 Hz where the differences become 
somewhat larger. 

Throughout our frequency band these RS modes have 
vertical seismic-gravity-gradient e-folding lengths Zggg = 
1/A: > £>! = 5 m. Thus, the upper parts of layer 2 
contribute significantly to the reduced transfer function 
/3', along with all of layer 1. 

For modes RS1-RS5, the gravity gradients are largely 
due to the S-waves' vertical surface motions, and corre- 
spondingly the reduced transfer functions have the fa- 
miliar range (3' ~ 0.6- 1.2 tha t we encountered for RS 
modes at Hanford (Sec. \ll B 3| ) and for SV-Up modes in 
a homogeneous half space (Fig. ^. 

By contrast, modes RS6-RS10 show a phenomenon not 
exhibited at Hanford: a broad dip in /3' to a value <C 1. 
This dip is caused by a significant excitation of P-waves in 
Livingston's 5-meter thick top layer: the vertical surface 
motions in the dip are largely due to the P-waves, and 
mass conservation guarantees that the gravity gradients 
they produce will be nearly cancelled by those from the 
subsurface, P-wave compressional source. These surface- 
layer excitations are not associated with any RP mode; 
as we shall see in the next subsection, the lowest RP 
mode at these frequencies has ch about twice as high as 
for these modes. It seems that the close proximity of the 
two very sharp geophysical discontinuities (the earth's 
surface and the sharp rise of cp caused by water) forces 
the modes' S-waves to generate a sizable component of 
P-waves even moderately far from P-wave resonance. No 
such phenomenon was observed in our Hanford 4-layer 
model. 



3. Mode RPl 



Figure ^ shows the dispersion relation for the lowest 
RP mode, RPl, at Livingston, along with the RF and 
lowest 10 RS modes. As noted earlier, RPl does not 
overlap the other modes [by contrast with Hanford (Figs. 
|and|)]. 



18 



1500 

o 

E 1250 




500 

o 

^ 250 







RP1 










RS10 










-.,RS7 
















RF _ 


5 10 


15 


20 


25 


3 



1.14 




5 10 15 20 25 30 

Frequency, f (Hz) 

FIG. 11. Upper panel: dispersion relation for mode RPl 
in the 4-layer Livingston model. Also shown for compari- 
son are the dispersion relations for the fundamental Rayleigh 
mode RF and the lowest few RS modes (c/. Fig. Lower 
panel: properties of mode RPl. 



pendix More specifically, the dispersion relation (with 
both ch and / real) was taken from the P-SV decoupling 
approximation, the P-waves with this ch and / were re- 
garded as impinging from layer 2 onto the top of layer 1 
at a glancing angle, and the reflected P- and SV- waves 
were regarded as propagating off to "infinity" (or, more 
realistically, dissipating before any return to the inter- 
face). This is the approximation that was so successful 
for the RS modes when combined with the correct 4-layer 
dispersion relation, but we don't have a good handle on 
its accuracy here, with the less reliable P-SV-decoupling 
dispersion relation. We are much more confident of our 
approximation for / > 22.8 Hz. There we used the exact 
two-layer equations (Appendix |c|), together with Lee's 
exact, complex dispersion relation ch (/)• 

These computations produced an anisotropy that 
peaks at f — 22.8 Hz where ch = cp2, with a peak 
value of ^ ^ 8 (Fig. ^l]). This is smaller than the 
peak anisotropics for mode RPl at Hanford (Fig. |^), but 
comparable to those for the higher-order RP modes that 
propagate nearly horizontally in the Hanford basalt (Fig. 

The reduced transfer function /?' lies in the same range, 
to 0.15, as for all the Hanford RP modes that we studied 
(Figs. and ^). This adds to our conviction that this low 
range of /?' is a general characteristic of RP modes. 



At frequencies / < 22.8 Hz, the RPl mode has hori- 
zontal speed Ch > cp2 — 1660 m/s and thus its P-waves 
can propagate in layers 1 and 2 (and also in layer 3 below 
11.3 Hz). In this region we have evaluated cuif) using 
the P-SV decoupHng approximation [Eq. (^)]. 

At frequencies / > 22.8 Hz, the horizontal speed is 
Ch < cp2, so the mode's P-waves are evanescent in layers 
2, 3, and 4. In this regime we have adopted an approxi- 
mation that is much more accurate than the decoupling 
one. We have idealized the material as two-layered: a 5 
m thick upper layer with the properties of layer 1 of Table 
IV, and below that a homogeneous half space with the 
properties of layer 2. For this layer-plus-half-space model 
we have used an analytic dispersion relation due to Lee 
po| (Appendix P). Because the mode's SV- waves can 
leak out of layer 1 into layer 2 (and then propagate away 
to "infinity" — or, more realistically, dissipate), Lee's 
dispersion relation predicts a complex frequency / if ch 
is chosen real, and a complex ch if / is chosen real. The 
predicted losses are small (quality factors Q decreasing 
from ~ 50 at / ~ 24 Hz to ~ 15 at 30 Hz). The real part 
of the dispersion relation is shown in the upper panel of 
Fig. [l|- 

The lower panel of Fig. |ll| shows the anisotropics 
and reduced transfer functions for this RPl mode. At 
/ < 22.8 Hz, where the P-waves are propagating nearly 
horizontally in layer 2, these properties were computed 
using the P-up approximation in the above two-layer 
(layer-plus-half-spacc) model; c/. the introduction to Ap- 



4. Summary of Livingston model results 

The most important of the above results are those for 
the reduced transfer functions /?' of the various modes at 
Livingston. They are summarized in Table | and their 



implications are discussed in the Introduction, Sec. ID 2 



V. CONCLUDING REMARKS 

A. Summary 

In this paper, we have used the theory of seismic sur- 
face waves to calculate the seismic gravity-gradient noise 
spectra that are to be expected at the Hanford, Wash- 
ington and Livingston, Louisiana LIGO sites. Our final 
noise strengths, as shown in Fig. 2, are close to Saul- 
son's previous rough estimate. At noisy times and near 
10 Hz, the seismic gravity-gradient noise is likely to be 
more serious than vibrational and thermal seismic noise 
in advanced interferometers. Unless means are found 
to combat gravity-gradient noise (see below for possible 
methods) , the hard- won gains in sensitivity due to R&D 
on vibration isolation and thermal noise may be com- 
promised by seismic gravity gradients, at least at noisy 
times. 
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B. Effects of topography and of LIGO construction 

In our analysis we have idealized the earth's surface 
near the LIGO test masses as perfectly planar and as 
undisturbed by LIGO construction. Irregularities in to- 
pography will significantly disturb the waves' propaga- 
tion and their vertical structure only if the surface height 
varies by amounts as large as ^ 2m/(//10Hz) = ('--^ 1/2 
the shortest vertical e-folding length Zp for RF waves), 
on horizontal lengthscales as short as ^ 8 m/ (//lOHz) ~ 
(~ 2 times the horizontal reduced wavelength 1/fc of 
those RF waves), within distances of the test masses 
~ 25m/(//10Hz) = (the horizontal wavelength of those 
RF waves), for frequencies ^3 — 30 Hz. (Of all the 
modes we have studied, the RF modes hug the surface 
most tightly and thus will be most influenced by the to- 
pography.) 

Variations on these scales were rare at the two LIGO 
sites before construction. However, the grading that 
made the arms flat produced topographic variations in 
the vicinity of some of the test masses that are marginally 
large enough to disturb the propagation. Examples are 
the long pits dug alongside the arms at Livingston to get 
material for building up the arms' heights, and excava- 
tion at Hanford to lower the arms below the level of the 
surrounding land near the southwest arm's midstation 
and the northwest arm's endstation. 

We speculate that these topographic modiflcations will 
alter the seismic gravity gradient noise by a few tens of 
percent, but probably not by as much as a factor 2. Fu- 
ture studies should examine this issue. 

The 1 m deep concrete foundations of the buildings 
that house the test masses will likely also influence the 
noise by a few tens of percent, particularly at '--^ 20 — 30 
Hz where the RF waves' vertical penetration is short. 
The foundation extends approximately 10 meters by 25 
meters at the interferometer's end stations (and also, in 
the case of Hanford, at the mid station). The foundation 
is approximately "X" shaped for the corner stations, with 
each arm of the "X" extending roughly 100 meters by 20 
meters |^l[. The sound speeds in the concrete will be 
a factor of several higher than the surrounding ground, 
so the foundations will form very sharp "geophysical" 
interfaces in the ground, causing diffraction of impinging 
waves and altering their vertical structure. Because the 
foundations are so shallow, we doubt that their net effect 
on the seismic gravity gradient noise can be as large as a 
factor 2, but future studies should examine it. 



C. Measurements that could firm up our 
understanding of seismic gravity gradients 

Our analysis is plagued by a large number of uncertain- 
ties regarding the true make-up of the ambient seismic 
background at the LIGO sites. We made extensive use 
of measurements of ground motion which functioned as 



constraints on what modes could be present. These mea- 
surements were helpful, but certain other measurements 
would be considerably more helpful. We suggest that, 
to the extent that resources permit, these measurements 
be included in future seismic surveys for gravitational- 
wave interferometer sites, including future surveys at the 
LIGO sites. 

First, we recommend careful measurements of the 
sound speeds and dynamical Poisson ratios of the ground 
as a function of depth, especially in the top few tens of 
meters and if possible down to the bedrock. At Hanford, 
we had reasonably complete data ||3^, thanks to earlier 
plans to build a nuclear power plant in the vicinity. As 
discussed in this paper, we encountered serious discrep- 
ancies between those old data and data from the LIGO 
geotechnical survey. At Livingston, we had no P-wave 
speed or Poisson ratio proflles, and the S-wave speed pro- 
files available only went down to a depth of 15 meters. 
As a result, we had to use a mixture of theory, profiles 
from other sites, and phenomenological fitting to obtain 
a plausible velocity profile. Velocity proflles are of cru- 
cial importance in determining how the various modes 
behave in the ground. 

Second, we recommend measurements that more 
nearly directly determine the modes that characterize 
the seismic motion. In this paper, as discussed above, 
we were able to put together very rough estimates of the 
modes that actually characterize the seismic background 
by using surface motion data as constraints, particularly 
anisotropy ratios measured at the sites, and by appealing 
to more detailed measurements at other sites. However, 
other techniques could provide much more useful and re- 
strictive constraints, thereby more sharply differentiating 
among the various modes. In particular: 



• Surface seismic arrays [£1|,^ allow one to measure 
the phase relationships of ground motion at appro- 
priately separated points, from which one can infer 
the excited modes' wave numbers k(f) and hori- 
zontal propagation speeds cnif)- 

• Borehole measurements |Q allow one to measure 
the phase correlation of motion at the surface and 
at some depth z underground, and the variation of 
amplitudes with depth, thereby introducing addi- 
tional constraints on the background. 

• Specialty seismic instruments called "dilatometers" 
[pS 29 1 measure directly the fractional density per- 
turbation Sp/p that is the subsurface source of seis- 
mic gravity gradients. Measurements down bore- 
holes with such devices could place further con- 
straints on the mode mixtures present, and could 
show how Sp/p varies with depth, at fixed fre- 
quency. When correlated with vertical surface seis- 
mic measurements, they could give information 
about the cancellation of gravity gradients from the 
surface and subsurface sources. 
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D. Mitigation of seismic gravity gradient noise 

Seismic gravity gradients are unlikely to be a major 
concern to LIGO detectors in the near future, since these 
detectors are only sensitive to frequencies / > 35 Hz. 
Eventually, however, LIGO experimenters may succeed 
in achieving extremely good vibration isolation and ther- 
mal noise control at frequencies / ^ 10 Hz. At this time, 
the detectors may well be plagued by seismic gravity- 
gradient noise, at least at noisy times; and there may be 
a strong need to try to mitigate it. 

We see two possibilities for modest amounts of mitiga- 
tion: (i) monitoring the noise and removing it from the 
LIGO data, and (ii) building moats to impede the prop- 
agation of RF-mode seismic waves into the vicinities of 
the test masses. 

Monitoring and correction: By using dedicated 3- 
dimcnsional arrays of vertical surface seismometers and 
borehole- mounted dilatometers in the vicinities of all test 
masses, one might be able to determine both the surface 
and subsurface components oi 5p/ p with sufficient spatial 
and temporal resolution for computing the seismic grav- 
ity gradient noise and then removing it from the data. 
This paper's insights into the modes at the two LIGO 
sites and the gravity gradients they produce may provide 
a foundation for future explorations of monitoring-and- 
correction strategies. 

Moats: By constructing a narrow, evacuated moat 
around each test mass, one might succeed in shielding 
out a significant portion of the RF waves that we sus- 
pect are the dominant source of quiet-time seismic grav- 
ity gradients. Since the RF mode contains substantial 
S-waves and they are the dominant contributors to the 
gravity-gradient noise, such moats may have to be at 
least as deep as the S-wav es' vertical e-folding length, 
Zs — 9-15m(10Hz//) [Eq. ( |2.5| ) modified for an increase 
in the RF speed ch due to stratification as shown in Figs. 
^ and ^. Since Zg — 2.5Z-p, moats of this depth would 
strongly shield out the RF mode's P-waves. 

The radius of the moats should be ^ A ^ 
20-35m(10Hz//). It is not clear to us whether such 
moats at Livingston would be effective if filled with wa- 
ter, or whether they would have to be kept pumped out. 
The water would shield out the RF mode's S-waves but 
transmit its P-waves. If, after transmission, the waves re- 
main mostly of P-type, then a significant reduction of (3' 
could result; but it is not at all obvious how much regen- 
eration of S-waves would occur in the moat-surrounded 
cavity. Detailed modeling would be required to sort out 
such issues. 

Although moats may be well-suited to reduce gravity 
gradients generated by the RF mode, they are probably 
not so well-suited to reduce gravity gradients generated 
by Rayleigh overtones. The overtones can be visualized 
as seismic waves that propagate by bouncing between 
layer interfaces and the earth's surface; they could prop- 
agate right under the moat and into the region under 



the test mass. Conceivably, they could even resonantly 
"ring" the earth under the mass, worsening the seismic 
gravity-gradient noise. 

If seismic gravity gradients become a problem in the fu- 
ture, ideas such as moats and monitoring-and-correcting 
will have to explored. 
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APPENDIX A: GENERAL EXPRESSION FOR 
REDUCED TRANSFER FUNCTION 



In this Appendix we derive Eqs. ( 1.17 )-( 1.24 ) for the 
reduced transfer function and anisotropy ratio of an ar- 
bitrary Rayleigh mode. In the text the mode is labeled 
J; in this Appendix we shall omit the subscript J. 

The mode has frequency /, angular frequency u = 
27r/, horizontal wave number k, horizontal phase speed 
Ch — oj/k, and horizontal propagation direction k. At 
the earth's surface its displacement vector is 



C(z = 0) = {(,Hk - ^ve,)e 



i{k-X — LOt) 



(Al) 



[Eq. ( 1.21 )] ; and on and beneath the surface it produces 
a fractional density perturbation 



6p 



= [Cy<5(z)-^7^(z)]e 



i{k-x—ujt) 



(A2) 



[Eq. ( [l.22| )]; here k = kk is the horizontal wave vector 
and S{z) is the Dirac delta function. 

Since the ambient seismic motions are horizontally 
isotropic, this mode is excited equally strongly for all hor- 
izontal directions k, and also for all wave numbers in some 
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(arbitrarily chosen) small band Afc around k — i.e., in the 
annulus C^k of width Afc in wave-vector space. Corre- 
spondingly (with an arbitrary choice for the strength of 
the excitation) , the net displacement along some horizon- 
tal direction n, in the frequency band A/ = c/fAfc/27r, 
is 



x{t) = n 



(A3) 



and the power of this random process X{t) in the fre- 
quency band A/ is 

X\f)Af = Yl \^^\"Ck ■ n? - \^H?^ , (A4) 



where N/^v — 1 the (normalization-dependent) to- 
tal number of allowed k values in the annulus Caa;, and 
the 1/2 comes from averaging (fc • n)'^ over the horizon- 
tal direction k. (Note: the overall normalization N^k 
of our procedure for going from the random process ex- 
pressed as a sum over directions to the processes's power 
will have no influence on our final answers for A and 
/3', since they are square roots of ratios of powers from 
which N^k drops out.) Similarly, the net displacement 
and power along the vertical direction are 



Z{t) = 



_^^^^tik-S-u^t) 



and 



Z'if)Af = J2\S.v\'^\^v\'K 



Ak 



(A5) 



(A6) 



The mode's anisotropy ratio, A~Z/Xis therefore 

A^V2\^v\mH\, (A7) 



cf. Eq. (1.23); and the direction-averaged power W^Af = 
{2X^Af + Z^Af)/3 is 



W^Af 



\^H\' + \iv\' 



N 



Ak 



(A8) 



By analogy with Eq. (A3), the isotropically excited 
mode produces a fractional perturbation in density on 
and beneath the earth's surface given by 



6p 
P 



3fJ 



i{k-x — uji) 



(A9) 



cf. Eq. ([A^). As an aid in computing the gravitational 
acceleration produced on one of the interferometer's test 
masses by these density perturbations, we place the ori- 
gin of coordinates (temporarily) on the earth's surface, 



immediately beneath the test mass. Then the location of 
the test mass is —Hez, where H is its height above the 
surface. We denote by m the unit vector along the laser 
beam that is monitoring the test mass's position. Then 
the gravitational acceleration along the m direction is 



am{t) = - 



3 , {x' ■rh)GSp{x',t) 



(AlO) 



Invoking Eq. ( A9) and introducing Cartesian coordinates 
{x',y',z') insi de th e sum with k along the x'-direction, 
we bring Eq. ( |A10| ) into the form 



k 



Gp 



{x'm, + y'my)e'''^ [iv5{z') -f TZjz')] 
\x'^ + y'2 + (z' + H)2]3/2 



dz'dx'dy' . 
(All) 



Integrating out the horizontal directions x' and y' from 
— oo to +00 at fixed 2', and integrating out the S func- 
tion, we obtain our final expression for the gravitational 
acceleration on the test mass 



E 

k 



2TTiGp{rh ■ fc)e"*"*e"'=^ 



(v+ I n{z')e-'''' dz' 



(A12) 



We next solve the pendular equation of motion for the 
displacement 6xj ■ rhj of the test mass in response to this 
gravitational acceleration (where the label = 1, 2, 3, or 
4 indicates which of the interferometer's four test masses 
we are discussing); the result is 



dxj ■ rrij 



2TriGp{k ■ rhj 

k 



'^^i{k-Xj—Ljt) ^ — kH. 



00 



X ( + niz')e 



dz' 



(A13) 



Here and t are the angular eigenfrequency and damp- 
ing time of the test mass's pendular motion. After com- 
pleting the calculation we have moved the origin of co- 
ordinates to the interferometer's beam splitter, thereby 
producing the term ik ■ Xj in the exponential, where Xj 
is the test mass's location; cf. Fig. |l2|. 

The interferometer's displacement signal x{t) = Lh{t) 
is its difference in arm lengths. 



x{t) 



4 

3 = 1 



(A14) 



We have chosen rhj to point away from the test mass's 
mirror on the first arm and toward the mirror on the 
second arm as shown in Fig. |l^. The seismic gravity- 
gradi ent n oise is obtained by inserting expression (A13) 
into ( A14) for each of the four test masses. 
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FIG. 12. The geometry of the interferometer. 



The contributions to this noise coming from the two 
end masses, j = 3 and 4, are not correlated with those 
coming from any other test mass in our 3-30 Hz fre- 
quency band, since 3 and 4 are each so far from the cor- 
ner and each other (4km A = 27r/fc). However, there 
is a significant correlation between the two corner test 
masses, 1 and 2. Taking account of this correlation, the 
i nterf erometer's displacement signal x{t) [Eqs. ( A13 ) and 
exhibits the following noise power in the frequency 



band A/: 



{2TTGpf 



E 

Afc 



/ ■R{z')e''''' dz' 



Jk, (A15) 



where 



ik-X2 |2 



-{k ■ rhzY + ■ '^T-iY 



(A16) 



into Eq. ( A17 ), setting k-rhx = coscj) and fc-TO2 = sin0 
and averaging the quantity inside the sum over k (i.e. 
over 0), we obtain 



(A18) 



where 



liv) 




cos (p sm (p cos y 



^/2 



(A19) 



[Eq. (|l.l9| )] . This function is graphed in Fig. ^ Inserting 
Eq. ( AISP into Eq. ( |A15[ ) and noting that ^^j. = 
= A^Afe is the number of allowed wave vectors in 
the annulus Caa;, we obtain our final expression for the 
interferometer's displacement noise power: 



{ATiGpf 



-2kn 



£,v+ I n{z')e-'''' dz' 



NAk 



. (A20) 



The transfer function T(/) for the seismic gravity- 
gradient noise is obtained by dividing the direction- 
averaged ground displacement noise power (A8) into the 
interferometer displacement noise power (A20) an d tak- 
ing the square root. The result is expression ( |l.8| ) with 
t he re duced transfer function f3 gi ven by (3 — ^T(3' [Eq. 
(ll.lTl) ], where T = e'^'^ [Eq. (|l.20| ) in which uj/chj = k] 
and 



/?'(/) 



3/2 



£,v+ I n{z')e-''''' dz' 





(A21) 



[Eq. (1.24)] 



APPENDIX B: FUNDAMENTAL RAYLEIGH 
MODE IN HOMOGENEOUS HALF SPACE 



Here we have broken up the sum over k into one over 
all directions k and one over its length k in the range 
Afc. Each of the last two terms in Jk (the uncorrelated 
contributions of masses 3 and 4) average to 1/2, and the 
first term can be rewritten in terms oi xi — X2'- 



Jk 



\k • jfiic 



-l-fc-mal^ + l 



(A17) 



By virtue of the geometry of the interferometer's corner 
test masses (Fig. |l^), xi — X2 = l{rhi + ?fi2)/\/2, where 
I is the separation between those masses. Inserting this 



In this Appendix we briefly review the theory of Rayleigh 
waves propagating in a homogeneous half space (i.e., a 
homogeneous, planar model of the earth), and then we 
derive the anisotropy ratio A and reduced transfer func- 
tion /?' for such waves. 

A homogeneous half space can support only the fun- 
damental Rayleigh mode, since the overtones all require 
inhomogeneities to confine them in the vicinity of the 
earth's surface. The theory of this mode is developed 
in a variety of standard texts [^-Q. According to 
that theory, the waves propagate with a horizontal speed 
ch which is slightly slower than the S-wave speed cs 
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(which in turn is slower than cp). The ratio Ch/cs 
is a function of the material's Poisson ratio varying 
from ch/cs — 0.904 for v — 0.16 (fused quartz) to 
ch /cs — 0.955 for v — 0.5 (fluids and other easily sheared 
materials). More generically, it is given by ch/cs — VQi 
where C, is the real root, in the range < C < Ij of the 
equation 



1 



c- 



= 



(Bl) 



The Rayleigh waves' horizontal wave number is fc = 
uj/ch, and their wavelength is A = 2'K/k. The P-wave 
of the fundamental Rayleigh mode decays with depth z 
as e"'^*'^, where the dimensionless ratio q of vertical e- 
folding rate to horizontal wave number is 



= v/1 - (cif/cp) 



(B2) 



Similarly, the SV-wave part decays with depth as e^*'^^, 
where the dimensionless ratio s of vertical e-folding rate 
to horizontal wave number is 



v/1 - [ch/cs? = 



(B3) 



More specifically, the mode's displacement eigenvector 
^ can be decomposed into a P-wave which is the gradient 
of a scalar potential plus an SV-wave which is the curl of 
a vector potential. We shall denote by -0 the complex am- 
plitude of the scalar potential. The normal components 
of elastodynamic stress produced by this wave must van- 
ish]^ at the earth's surface. Upon imposing these bound- 
ary conditions, a standard calculation ||3l| , ^^ gives the 
following expression for the displacement vector: 



-qkip I e 



,-qkz 



1 + s^ 



(B4) 

Here, is the unit vector pointing in the z-direction, 
which we take to be down, t is time, x denotes hori- 
zontal location, and k = kk is the mode's horizontal 
wave vector. By comparing this displacement vector with 
Eq. (1.21), we read off the following expressions for the 
horizontal and vertical displacement amplitudes at the 
earth's surface, z = 0: 



The wave displacement (B4) produces a fractional per- 
turbation 5p/p of the earth's density beneath the surface 
given by 



6p{z > 0) 



where 



n{z) = (1 - q^)k'''tPe 



(B6) 



(B7) 



Inserting Eqs. (^) into Eq. (|]2|), we obtain the 
anisotropy ratio for the RF mode of a homogeneous half 
space, 



q{l - s^) 



2qs 



(B8) 



and inserting ( |B5| ) and (B7) into ( 1.24 ) and integrating, 
we obtain the mode's reduced transfer function 



3(1 



2qy 



2(1 + s^)[il + S^)il + q^) - Aqs] 



(B9) 



1 



2qs 



1 + S2 



(B5) 



^More accurately, they must be continuous with the stress 
produced by the earth's atmosphere, which we approximate 
as vacuum. 
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APPENDIX C: MULTILAYER MODEL 



In this Appendix we derive the equations governing Rayleigh overtones and the reduced transfer function in a multi- 
layer model of geophysical strata. 

1. Model and notation 

Our model consists of N homogeneous layers labeled by the index n = 1,2,3, ... ,N. Layer 1 is at the surface, 
layer A'' is a homogeneous half space at the bottom, and the interfaces between layers are horizontal. The Rayleigh 
modes propagate as decoupled planar SV- and P-waves in each layer; they are coupled at the interfaces by continuous- 
displacement and continuous-normal-stress boundary conditions. 

We have already introduced much of our notation in the body of the paper; to make this Appendix self-contained, 
we reiterate some of it here, along with some new notation: 

u! = 27r/: Angular frequency of waves. 

k = kk: Horizontal wave vector, with k its magnitude and k the unit vector in its direction. 

Ch = u)/k: Horizontal phase velocity of waves. 

Cz'. Downward pointing unit vector. 

Dn'. Thickness of layer n. 

Zn- Depth below the top of layer n. 

^n- Displacement vector for waves in layer n. 

Kn- Bulk modulus in layer n. 

Hn- Shear modulus in layer n. 

Pn'. Density in layer n. 

cpn- Speed of propagation of P-waves in layer n. 
csn- Speed of propagation of S- waves in layer n. 

apn- Angle to vertical of P-wave propagation direction (between and 7r/2 if real, by convention). If P-waves are 
evanescent in the layer, ap„ will be complex. 

asn- Angle to vertical of SV-wave propagation vector (between and 7r/2 if real, by convention). If SV- waves are 

evanescent in the layer, agn will be complex. 

Vn- Complex amplitude of upgoing P-waves at the top of layer n. 
V'n- Complex amplitude of downgoing P-waves at the top of layer n. 
<S„: Complex amplitude of upgoing SV- waves at the top of layer n. 
S'n- Complex amplitude of downgoing SV- waves at the top of layer n. 

In accord with this notation, the displacement vector in layer n has the following form: 

+ (5;e*'="""°*"^'' - 5„e-^'=^"'^°*"^") cosas„ k - (^^e*'^^" ^ ^^g-ife^„ cot«s„) ^y^asn 4 

"(CI) 

Since the waves are generated at the Earth's surface, the upward propagating waves are absent in the lowermost layer: 

Pjv = 0, Sn=Q. (C2) 
Consequently, the waves have 4 A'' — 2 complex amplitudes. 
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2. Equations for the dispersion relation, the propagation angles, and the amplitudes 



Once one has specified the Rayleigh mode of interest, its horizontal propagation direction fc, and one of its am- 
phtudes, say Vi, then all its other properties are uniquely determined as a function of frequency. To evaluate its 
properties one first computes its horizontal dispersion relation uj{k) [or equivalently cnif)] by a procedure to be 
outlined below. Then one computes all the waves' propagation angles by imposing Snell's law (i.e., by demanding 
that all components of the wave propagate with the same horizontal speed ch)'- 

CH ■ (C3) 



sm a pn sm asn 

At the Earth's surface, the (primed) amplitudes of the reflected waves are related to the (unprimed) amplitudes of 
the incident waves by the following two standard equations ||3l|^^3^ : 

2 sinasi cosa_pi('P{ — Vi) + cos2a5i(5( + Si) = 
sinapi cos2asi(P( + "Pi) — sinagi sin2asi(iS{ — Si) 

= . (C4) 

These equations can be derived by setting the vertical- vertical and vertical-horizontal components of the stress to zero 
at the Earth's surface, and by expressing the ratio of bulk to shear modulus in terms of the propagation angles: 



(C5) 



Kn _ cpn"^ ^ 4 _ sin^ apn _ 4 
csn^ 3 sin^ asn 3 

The junction conditions at the interface between layer n and layer n -I- 1 take the following form p3| , ^ : 

{Vy^^- cotap„ _^ -p^e-'^D^ cotap„^ siu ap„ + {S'^'^^^ - SnC-'^^^ cotas„) cOS as„ 

= (^n+i + Pn+i) sinap„+i -f (5^+1 - 5„+i) cosa5„+i , (C6a) 

^p/^gifcD„cotap„ _ p^g-»/cD„cotap„^ COS ap„ - (^^g*^- ™t as„ + 5^^g-*r.„ cot as„) gin as„ 

= ^ Pn+i) cosap„+i - {S'n+i + S^+i) sina5„+i , (C6b) 

^in [(1 - COi'aSn) (75;e'fcl3,.cotap„ ^p^^g-,fcD„cotap„-) ^j^^ ^^^^ ^ 3 (5;e'fc^" -t as„ _ ^^g-.fcl?,. cot as„ ) cOS agn] 
= ^„+i [(1 - COt^ asn+l) i'P'n+l + 'Pn+l) siu ap„+i -I- 2 (5^+1 - Sn+l) COS as„+i] , (C6c) 

fin[2[Pne " ^" -VnG ^" j COS ap„ - ( 1 - cot asn) [S,^e " '"^+Sne " """jsmasnj 

= ^„+i [2 {V'n+i - Vn+l) cosapn+i - (1 - cot^ aSn+l) {S'n+i + Sn+l) sinQ;5„+i] . (C6d) 



Equation ( |C6a| ) is continuity of the horizontal displacement, ( C6b|) is continuity of the vertical displacement, ( |C6c[ ) is 



continuity of the vertical- vertical component of the stress, and (C6d) is continuity of the vertical- horizontal component 
of the stress. 



Equations (C4) and (p6a|)-(C6d) are AN — 2 homogeneous linear equations for AN — 3 independent ratios of 



amplitudes, and for the horizontal dispersion relation uj{k) [or equivalently cnif)]- It is convenient to evaluate the 



dispe rsion relation by setting to zero the determinant of the coefficients of the amplitudes in Eqs. ( p4| ) and (C6a) 
( |C6d| ). The remaining AN — 3 amplitudes can then be computed in terms of Vi using any AN — 3 of these equations 



This was the procedure used to derive the 4-layer results quoted in the text. Once the dispersion relation and the 
amplitudes have been evaluated as functions of frequency, the anisotropy ratio and reduced transfer function can be 
computed using the equations derived in the following subsection. 

3. Anisotropy ratio, and reduced transfer function 



By comparing Eq. (1.21) with the displacement eigenfunction (pl[) for layer n = 1, we read off the horizontal and 



vertical displacement amplitudes at the earth's surface: 
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= {V-^ +'Pi)sinapi + {S[ -Sijcosasi , 



The wave displacement (CI) produces a fractional density perturbation 5pn/ Pn 
layer n, with amplitude given by 



sm a pr 



I ikZji cot Qpn 



-ikzn cot apTi ^ 



By inserting Eqs. (p7|) and (jCq) into Eq. (1.23), we obtain the anisotropy ratio 



A = V2 



{V'l - Vi) COS api - {S[ + Si) sin 0151 



{V-^ + Pi) sinapi + {S[ - 5i) cosagi 



By inserting Eqs. (CTp, (C8), (Cg), and the relation 



n-l 



(C7) 
(C8) 

-V • in = 7^„(z„)e'('■■"~"*^ in 

(C9) 

(CIO) 
(Cll) 



into Eq. (1.24), integrating, and summing over all four layers, we obtain the reduced transfer function 

A^(/) 



/3'(/) = 



where 



3 

N 



{Vi — v'l) cosap + {Si + S[) sin as 



n—l 



-[kD„{l-i cot ap^)] 



(C12a) 



(C12b) 



(T'^ + T'l ) sin api + (5( - 5i ) cos asi 



{V'l - Vi) cos api - {S'l + Si) sin asi 



(C12c) 



In Eq. (C12b) for Af{f), we have inserted the factor Pn/pi to allow for the possibility (ignored in the text) that the 
different layers have different densities. 
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APPENDIX D: LEE'S DISPERSION RELATION 
FOR 2-LAYER MODEL 



When there are only two layers, a top layer with thickness 
D and a bottom layer with infinite thickness, the disper- 
sion relation u}{k) [or equivalently c/f(/)] of the multi- 
layer model (Appendix y) can be brought into an ex- 
plicit form that permits rapid numerical solutions. This 
form was derived by Lee pO[ by manipulating the 6x6 
dete rmina nt of t he co efficients of the amplitudes in Eqs. 
( p^ and ( |C6aD -~( |C6d| ). The standard textbook by Erin- 
gen and §uhubi |34[| presents and discusses Lee's disper- 
sion relation [pages 547-550; note that on the first line 
of their Eq. (7.7.44) V2 should be vi\. The dispersion 
relation consists of the following prescription: 
The unknown to be solved for is 



C = {CH/Cs2f 



(Dl) 



At low propagation speeds ch (high frequencies) the SV- 
waves in layer 1 will typically propagate rather than de- 
cay, with vertical wave number divided by horizontal 
wave number given by 



= ^/C{cs2/csiY - 1 = V [ch/csiY ~ 1 

— cot asi , (D2a) 

while the other waves will typically be evanescent with 
ratios of e-folding rate to horizontal wave number given 
by 

qi = Vl-C(cs2/cpi)2 = ^\-{cHlcpif, (D2b) 

92 = V1-C(CS2/CP2)2 = ^1-(ch/cp2)', (D2c) 

S2 = = ^Jl-{cH/cs2)\ (D2d) 

Regardless of the magnitude of ch and thence regard- 
less of whether these quantities are real or imaginary, we 
regard them all as functions of ch given by the above 
expressions. 

We define two quantities 



Q = A*2/M1 , R = Pl/P2 



(D3) 



that appear in what follows. In terms of C, Q, and R, we 
define 



X = QC-2(Q-1), 
Y = QRC + 2{Q-l) , 
Z = Q(l-i?)C-2(Q-l) , 
W = 2{Q-1) . 



(D4) 
(D5) 
(D6) 
(D7) 



In this dispersion relation and only here A, Y, Z, W rep- 
resent these functions instead of representing earth dis- 
placements. In terms of the above quantities we define 



a = (1 - ^i') 

+ 2(71 

6 = (1 - ai') 

+ 2(71 

77i = (l-^i') 
+ 2(?i - 

+ 2gi 



Xcosh{kqiD) + —Ysmh{kqiD) 
9i 



q2W sm{kaiD) Z cos{kaiD) 

0-1 



S2W cosh{kqiD) H Z smh{kqiD) 

I1 



Xsm{kaiD) - —Y cos{kaiD) 
CTi 



q2W cos{kaiD) H Zsm{kaiD) 



(D8) 



(D9) 



Asinh(fcgiD) - — F cosh(/cgii:i) 
9i 



Xcos{k(JiD) + —Ysm{kaiD) 
^1 



S2W^sinh(fc9ii:)) Z cosh{kqiD) 
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(DIO) 



(Dll) 



In terms of these four quantities. Lee's dispersion relation 
takes the form 



F(C, kD) = ei?72 - 6m = 



(D12) 



In the language of Lee's dispersion relation, finding 
multiple Rayleigh modes is a m atter of finding multiple 
curves C{kD) that satisfy ( D12 ). Each such C{kD) can be 
translated into a corresponding cni f), since ch — \JXcs2 
and / — cukjl'K. Overtone modes undergo a transition 
in layer 2 from propagating and lossy (so that seismic 
wave energy is lost from layer 1 into layer 2), to evanes- 
cent and confined (so the waves are restricted to the 
vicinity of the top layer)] at speed cnif) — cs2, which 
is equivalent to C, — 1. Thus, to produce dispersion re- 
lations for overtone modes, one can look for solutions to 
( D12 ) in the vicinity of ^ = 1, and then, depending on 
whether one wants confined modes or lossy modes, trace 
them from C = 1 to higher frequencies and lower horizon- 
tal speeds, or to lower frequencies and higher horizontal 
speeds. 



In Sec. IV C 3 we use Lee's dispersion relation to study 



the RPl mode at Livingston in the lossy regime. 



APPENDIX E: P-UP AND SV-UP MODES 



In the text we encounter situations in which one can 
approximate an overtone mode as P- or SV-waves that 
propagate upward through a homogeneous half space un- 
til they encounter the Earth's surface or one or more lay- 
ers near the surface, and then (exciting the layers) reflect 
back downward with accompanying production of the 
other type of wave. Such "P-up" and "SV-up" modes can 
be described by the multilayer equations of Appendix |^, 
with the up (unprimed) amplitudes in the bottom layer 
(homogeneous half space), 5, set to {Vb 7^ 0,5b = 0} for 
P-up modes, and {Vb — Q,Sb^ 0} for SV-up modes. 



28 



We can derive simple formulas for the anisotropy ratio 
A and reduced transfer function f}' of such modes for the 
case of no surface layers (a pure homogeneous half space) : 



1. P-up modes in a homogeneous half space 



as = a,rcsm{cs / ch) , ap = arcsin(cp/cH) 



The displacement function is given by Eq. (CI) with 
the subscript n's deleted since there is only one layer. 
The primed (down) amplitudes are given in terms of the 
unprimc d (u p) amplitude V by the surface junction con- 
ditions (p4); in particular 



S' = 



4 cos ap sin as cos as — sin ap cos^ 2as , 

o I 

4 cos ap sin as cos as + sin ap cos^ 2as 

4 sin ap cos ap sin as cos 2as 

q 1 

4 cos ap sin as cos as + sin ap cos^ 2as 



(Ela) 
(Elb) 



Inserting these into Eq. (CIO) we obtain the following 
anisotropy ratio: 



A ^ \/2cot2as 



where, by Snell's law [Eq. (3.1)], 



as = arcsin(cs/c_f/) 



(E2) 



(E3) 



Inserting e xpre ssions (El) into the one-layer version of 
equations ( |C12 ), we obtain the following reduced transfer 
function: 



P' = \/6 sin^ as 



(E4) 



The anisotropics and reduced transfer functions of Eqs. 
(E2) and (E4) are shown graphically in Figs. || and ^ 
for cp/cs — 2 (approximately appropriate to the surface 
materials at Livingston and Hanford). 



2. SV-up modes in a homogeneous half space 

For SV-up modes, as for P-up modes, the displace- 
ment function is given by Eq. (CI) with the subscript 
n's deleted. The primed (down) amplitudes are given in 
terms of the unprimed (up) amplitude S by the surface 
junction conditions (C4); in particular 



v 



sin as sin 4as 



4 cos ap sin as cos as + sin ap cos^ 2as 
^, 4 cos a p sin'^ as cos as — sin a p cos^ 2as ^ 
4 cos a p sin'^ as cos as -I- sin a p cos^ 2as 



S, (E5a) 



(E5b) 



Inserting these into Eq. ( plO ), we obtain the following 
anisotropy ratio 



^ = 2V2 
where, by Snell's law. 



cot ap 



cot^ as 



(E6) 



(E7) 



Inserting expressions ( |E5| ) into the one-layer version of 
equations ( pl2 ), we obtain the following reduced transfer 
function: 

, \/6 sin^ as 1 1 — 2i cot ap sin^ as sec 2as I , , 

P = ; = ■ (-CjS) 

1 + (2 I cot ap I sin^ as sec 2as)^ 



T he anisotropics and reduced transfer functions of Eqs. 
(E6) and (|e|) are shown in Figs. || and^, for cp/cs — 2. 
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